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Abstract The purpose of this review article is to give an up to date account of the 
theory and application of scale functions for spectrally negative Levy processes. 
Our review also includes the first extensive overview of how to work numerically 
with scale functions. Aside from being well acquainted with the general theory of 
probability, the reader is assumed to have some elementary knowledge of Levy pro- 
cesses, in particular a reasonable understanding of the Levy-Khintchine formula and 
its relationship to the Levy-Ito decomposition. We shall also touch on more general 
topics such as excursion theory and semi-martingale calculus. However, wherever 
possible, we shall try to focus on key ideas taking a selective stance on the technical 
details. For the reader who is less familiar with some of the mathematical theo- 
ries and techniques which are used at various points in this review, we note that 
all the necessary technical background can be found in the following texts on Levy 
processes; Bertoin (1996), Sato (1999), Applebaum (2004), Kyprianou (2006) and 
Doney (2007). 
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Chapter 1 
Motivation 



1.1 Spectrally negative Levy processes 

Let us begin with a brief overview of what is to be understood by a spectrally neg- 
ative Levy process. Suppose that (£2,&,F,P) is a filtered probability space with 
filtration F = : t > 0} which is naturally enlarged (cf. Definition 1.3.38 of [17]). 
A Levy process on this space is a strong Markov, F-adapted process X = {X, : t > 0} 
with cadlag paths having the properties that P(Xq = 0) = 1 and for each < s < t, 
the increment X, — X s is independent of and has the same distribution as X t - S . In 
this sense, it is said that a Levy process has stationary independent increments. 

On account of the fact that the process has stationary and independent incre- 
ments, it is not too difficult to show that there is a function f:M^C, such that 

E{e wx >)=e-'^ e \ (>0,eeR, 

where E denotes expectation with respect to P. The Levy-Khintchine formula gives 
the general form of the function f. That is, 

'F(0)=i J U0 + ^0 2 + / (l-e lto + i0xl w<1 )n(dx), (1.1) 

for every 9 <G R, where fi € t, ff e 1 and fl is a measure on M\{0} such that 
/(l Ax 2 )fl(dx) < oo. Often we wish to specify the law of X when issued from 
x € K. In that case we write P x , still reserving the notation P for the special case Po- 
Note in particular that X under P x has the same law as x + X under P. The notation 
E x will be used for the obvious expectation operator. 

We say that X is spectrally negative if the measure fl is carried by (— °°,0), 
that is fl(0,°°) = 0. We exclude from the discussion however the case of monotone 
paths. Under the present circumstances that means we are excluding the case of a 
descending subordinator and the case of a pure positive linear drift. Included in the 
discussion however are descending subordinators plus a positive linear drift and a 
Brownian motion with drift. Also included are processes such as asymmetric a- 
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stable processes for a G (1,2) which have unbounded variation and zero quadratic 
variation. These processes are by no means representative of the true variety of 
processes that populate the class of spectrally negative Levy processes. Indeed in 
the forthcoming text we shall see many different explicit examples of spectrally 
negative Levy processes. Moreover, by adding independent copies of any of the 
aforementioned, and/or other, spectrally negative processes together, the resulting 
process remains within the class of spectrally negative Levy process. Notationally 
we say that a Levy process X is spectrally positive when —X is spectrally negative. 

Thanks to the fact that there are no positive jumps, it is possible to talk of the 
Laplace exponent y/(A) for a spectrally negative Levy process, defined by 

E(e**)= e *W, (1.2) 

for at least X > 0. In other words, taking into account the analytical extension of the 
characteristic exponent, we have y/(A) = — V(— iA). In particular 

W {X) = -piX + — X 2 + / (e ljc -l-Axl w<1 )n(dx). (1.3) 

l J (-oo,0) 

Further, using Holder's inequality and the fact that J(_„ j(l Ax 2 )n(dx) < °°, it is 
easy to check that y/ is strictly convex and tends to infinity as X tends to infinity. 
This allows us to define for q > 0, 

&(q) = sup{X > : y/(A) = q}, 

the largest root of the equation y/(A) = q. Note that there exist two roots when 
q = (there is always a root at zero on account of the fact that y/(0) = 0) and 
precisely one root when q > 0. Further, from a straightforward differentiation of 
(1.2), we can identify y/'(0 + ) = E(X\) € [— °°,°°) which determines the long term 
behaviour of the process. Indeed when ±y/(0+) > we have lim^X t = ±°°, that is 
to say that the process drifts towards ±<», and when y/(0+) = then limsup^X, = 
— liminf ( |ooX r = oo, in other words the process oscillates. 



1.2 Scale functions and applied probability 

Let us now turn our attention to the definition of scale functions and motivate the 
need for a theory thereof. 

Definition 1.1. For a given spectrally negative Levy process, X, with Laplace ex- 
ponent y/, we define a family of functions indexed by q > 0, : R — > [0,°°), as 
follows. For each given q > 0, we have W^ q \x) = when x < 0, and otherwise 
on [0,oo), ffW is the unique right continuous function whose Laplace transform is 
given by 
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for j3 > <P(q). For convenience we shall always write W in place of W^°\ Typically 
we shall refer the functions as q-scale functions, but we shall also refer to W 
as just the scale function. 

From the above definition alone, it is not clear why one would be interested in such 
functions. Later on we shall see that scale functions appear in the vast majority of 
known identities concerning boundary crossing problems and related path decompo- 
sitions. This in turn has consequences for their use in a number of classical applied 
probability models which rely heavily on such identities. To give but one immediate 
example of a boundary crossing identity which necessitates Definition 1.1, consider 
the following result, the so-called two-sided exit problem, which has a long history; 
see for example [107, 101, 89, 15, 16]. 

Theorem 1.2. Define 

T+ = inf{f > : X t > a} and X = inf{f > : X, < 0}. 
For all q>0, and x < a, 

In fact, it is through this identity that the 'scale function' gets its name. Possibly 
the first reference to this terminology is found in Bertoin [13]. Indeed (1.5) has 
some mathematical similarities with an analogous identity which holds for a large 
class of one-dimensional diffusions and which involves so-called scale functions 
for diffusions; see for example Proposition VII, 3.2 of Revuz and Yor [87]. In older 
Soviet-Ukranian literature is simply referred to as a resolvent, see for example 
[57]. 

In the forthcoming chapters we shall explore in more detail the formal analytical 
properties of scale functions as well as look at explicit examples and methodology 
for working numerically with scale functions. However let us first conclude this 
motivational chapter by giving some definitive examples of how scale functions 
make a non-trivial contribution to a variety of applied probability models. 

Optimal stopping: Suppose that x, q > and consider the optimal stopping problem 

v(x) = sup£(e-" T+ ^ Vx )), 

T 

where the supremum is taken over all, almost surely finite stopping times with re- 
spect to X andX, = sup s<r X v . The solution to this problem entails finding an expres- 
sion for the value function v(x) and, if the supremum is attained by some stopping 
time T*, describing this optimal stopping time quantitatively. 

This particular optimal stopping problem was conceived in connection with the 
pricing and hedging of certain exotic financial derivatives known as Russian options 



6 



Alexey Kuznetsov, Andreas E. Kyprianou and Victor Rivero 



by Shepp and Shiryaev [94, 95]. Originally formulated in the Black-Scholes setting, 
i.e. the case that X is a linear Brownian motion, it is natural to look at the solution 
of this problem in the case that X is replaced by a Levy process. This is of particu- 
lar pertinence on account of the fact that modern theories of mathematical finance 
accommodate for, and even prefer, such a setting. As a first step in this direction, 
Avram et al. [9] considered the case that X is a spectrally negative Levy process. 
These authors found that scale functions played a very natural role in describing the 
solution (v, T*). Indeed, it was shown under very mild additional conditions that 

v(x) = e x "w®(y)dyj 

where, thanks to properties of scale functions, it is known that the constant x* e 
[0,°°) necessarily satisfies 

x* =inf{x>0: l+q f* W M (y)dy ~ qW (q) (x) < 0}. 
Jo 

In particular an optimal stopping time can be taken as 

T* = inf{f > : (iVf,) -X, > x*}, 

which agrees with the original findings of Shepp and Shiryaev [94, 95] in the diffu- 
sive case. 

Ruin theory: One arm of actuarial mathematics concerns itself with modelling the 
event of ruin of an insurance firm. The surplus wealth of the insurance firm is often 
modelled as the difference of a linear growth, representing the deterministic col- 
lection of premiums, and a compound Poisson subordinator, representing a random 
sequence of i.i.d. claims. It is moreover quite often assumed that the claim sizes are 
exponentially distributed. Referred to as Cramer-Lundberg processes, such models 
for the surplus belong to the class of spectrally negative Levy processes. 

Classical ruin theory concerns the distribution of the, obviously named, time to 
ruin, 

T " :=inf{?>0:X f <0}, (1.6) 

where X is the Cramer-Lundberg process. Recent work has focused on more com- 
plex additional distributional features of ruin such as the overshoot and undershoot 
of the surplus at ruin. These quantities correspond to deficit at ruin and the wealth 
prior to ruin respectively. A large body of literature exists in this direction with many 
contributions citing as a point of reference the paper Gerber and Shiu [44]. For this 
reason, the study of the joint distribution of the ruin time as well as the overshoot 
and undershoot at ruin is often referred to as Gerber-Shiu theory. 

It turns out that scale functions again provide a natural tool with which one may 
give a complete characterisation of the ruin problem in the spirit of Gerber-Shiu 
theory. Moreover, this can be done for a general spectrally negative Levy process 
rather than just the special case of a Cramer-Lundberg process. Having said that, 
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scale functions already serve their purpose in the classical Cramer-Lundberg set- 
ting. A good case in point is the following result (cf. [18]) which even goes a little 
further than what classical Gerber-Shiu theory demands in that it also incorporates 
distributional information on the infimum of the surplus prior to ruin. 

Theorem 1.3. Suppose that X is a spectrally negative Levy processes. For t > 
define X, = inf s < t X s . Let q,x > 0, «, v,y > 0, then 



E x (e-«*o ■ -x z - g du, X x -_ e dv, X z -_ e dy) 

= l{o< > .<vA,, M >o } e-* ( '' )(v - y) {wW'(x-y)-0( 9 )wW(x-y)}n(-d«-v)dydv 



Note that the second term on the right hand side corresponds to the event that ruin 
occurs by creeping. That is to say the process X enters (— °°, 0) continuously for the 
first time. 

The use of scale functions has somewhat changed the landscape of ruin theory, 
offering access to a number of general results. See for example [86, 6, 72, 66], as 
well as the work mentioned in the next example. 

Optimal control: Suppose, as usual, that X is a spectrally negative Levy process 
and consider the following optimal control problem 



where x, q > 0, a n = inf{f > : X t - Lf < 0} and L 71 = {Lf : t > 0} is the control 
process associated with the strategy %. The supremum is taken over all strategies 
% such that L % is a non-decreasing, left-continuous adapted process which starts 
at zero and which does not allow the controlled process X — L K to enter (— oo,0) 
as a consequence of one of its jumps. The solution to this control problem entails 
finding an expression for u(x) and, if the supremum is attained by some strategy n* , 
describing this optimal strategy. 

This exemplary control problem originates from the work of de Finetti [40], again 
in the context of actuarial science. In the actuarial setting, X plays the role of a 
surplus process as usual and L K is to be understood as a dividend payment. The 
objective then becomes to optimize the mean net present value of dividends paid 
until ruin of the aggregate surplus process occurs. 

A very elegant solution to this problem has been obtained by Loeffen [75], 
following the work of Avram et al. [8], by imposing some additional assump- 
tions on the underlying Levy process, X. His solution is intricately connected to 
certain analytical properties of the associated scale function W^ q \ By requiring 
that, for x > 0, — JT(— x) has a completely monotone density (recall that 17 
is the Levy measure of X), it turns out that for each q > 0, the associated g-scale 



+ y (w^'(x) - 5 ( o,o,o) (d«,dv,dy) 



where 5( ,o.o)(dM,dv,dy) is the Dirac delta measure. 
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function, W^ q \ has a first derivative which is strictly convex on (0,°°). Moreover, 
the point a* := inf{a > : W {q) '{a) < W {q) '{x) for all x > 0} provides a threshold 
from which an optimal strategy to (1.7) can be characterized; namely the strategy 
L* = (a*VX t ) -a*,t>0. The aggregate process, 



thus has the dynamics of the underlying Levy process reflected at the barrier a*. The 
value function of this strategy is given by 



The methodology that lead to the above solution has also been shown to have appli- 
cability in other related control problems. See for example [76, 77, 71, 78, 68]. 

Queuing theory: The classical M/G/l queue is described as follows. Customers 
arrive at a service desk according to a Poisson process with rate X > and join a 
queue. Customers have service times that are independent and identically distributed 
with common distribution F concentrated on (0,°o). Once served, they leave the 
queue. 

(x) 

The workload, Y t , at each time t > 0, is defined to be the time it will take a 
customer who joins the back of the queue at that moment to reach the service desk 

(x) 

when the workload at time zero is equal to x > 0. That is to say, Y, is the amount of 
processing time remaining in the queue at time t when there was an initial workload 
x. If {X t : t > 0} is the difference of a pure linear drift of unit rate and a compound 
Poisson process, with rate X and jump distribution F, then it is not difficult to show 
that when the current workload at time t = is x > 0, then 



Said another way, the workload has the dynamics of a spectrally positive Levy pro- 
cess reflected at the origin. 

Not surprisingly, many questions concerning the workload of the classical M/G/l 
queue therefore boil down to issues regarding fluctuation theory of spectrally neg- 
ative Levy processes. A simple example concerns the issue of buffer overflow ad- 
justment as follows. Suppose that there is a limited workload capacity in the queue, 
say B > 0. When the workload exceeds the buffer level B, the excess over level B 
is transferred elsewhere so that the adjusted workload process, say {y f ' x B ' : t > 0}, 
never exceeds the amount B. We are interested in the busy period for this M/G/l 
queue with buffer overflow adjustment. That is to say, 



a*+X t -(a*VX t ), t >0, 




(1.8) 



Y t [x > = (xVX t )-X t ,t>0. 




inf{f > : Y t (x ' B) = 0}. 
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Once again, scale functions prove to be an instrumental tool. For the case at hand, it 
is known that for all q > 0, 

l+qf B WM(y)dy ' 

See Pistorius [82] and Dube et al. [36] for more computations and applications in 
this vein. 



Continuous state branching processes: A [0,°°] -valued strong Markov process 
Y = {Y t : t > 0} with probabilities {P x : x > 0} is called a continuous-state branching 
process if it has paths that are right continuous with left limits and its law observes 
the following property. For any t > and y\,y 2 > 0, Y, under P yi + y2 is equal in law 

to the independent sum f/ 1 ' + f/ 2 ^ where the distribution of f/'' is equal to that of 
Y t under P yi fori— 1,2. 

It turns out that there is a one-to-one mapping between continuous state branch- 
ing processes whose paths are not monotone and spectrally positive Levy processes. 
Indeed, every continuous-state branching process F may be written in the form 

Y t =X diATo ,t>0, (1.9) 
where X is a spectrally positive Levy process, = inf{f > : X t < 0} and 

6>,=inf{s>0: f^>t}. 
Jo x u 

Conversely, for any given spectrally positive Levy process X, the transformation on 
the right hand side of (1.9) defines a continuous-state branching process. (In fact 
the same bijection characterises all the continuous-state branching processes with 
monotone non-decreasing paths when X is replaced by a subordinator). 

A classic result due to Bingham [19] gives (under very mild conditions) the law of 
the maximum of the continuous-state branching process with non-monotone paths 
as follows. For all x > y > 0, 



W(x-y) 
W(x) ' 



P,,(supFs <x) = 



where W is the scale function associated with the underlying spectrally negative 
Levy process in the representation (1.9). See Caballero et al. [23] and Kyprianou 
and Pardo [69] for further computations in this spirit and Lambert [74] for further 
applications in population biology. 



Fragmentation processes: A homogenous mass fragmentation process is a Markov 
process X := {X(f) : t > 0}, where X(f) = (X\ (t),X 2 (t), ■ ■ ■ ), that takes values in 
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y := \ s = (ji,,S2,- ••) :si >*2> - > 



o,£j/<i| 



Moreover X possesses the fragmentation property as follows. Suppose that for 
each s £ P s denotes the law of X with X(0) = s. Given, for t > 0, that 
X(f) = (5i,52,- • •) € y, then for all u > 0, X(t + u) has the same law as the variable 
obtained by ranking in decreasing order the elements contained in the sequences 
X^\u) ,X^ 2 \u) , where the latter are independent, random mass partitions with 
values in having the same distribution as X{u) under P( ill o,o,-)> P(s 2 ,0,0,-)' ' ' ' 
respectively. The process X is homogenous in the sense that, for every r > 0, the 
law of {rX(t) : t > 0} under P (10 ...) is P( rj0 ,o,-)- 

Similarly to branching processes, fragmentation processes have embedded ge- 
nealogies in the sense that a block at time t is a descendent from a block at time 
s < t if it has fragmented from it. It is possible to formulate the notion of the evo- 
lution of such a genealogical line of descent chosen 'uniformly at random'. To 
avoid a long, detailed exposition, we refrain from providing the technical details 
here, and mention instead that if {X*(t) : t > 0} is the sequence of embedded frag- 
ments along the aforesaid uniformly chosen genealogical line of descent, then it 
turns out that {— logX*(t) : t > 0} is a subordinator. We suppose that, for 9 > 0, 
0(0) = E(i oo j ...)(^*(f) 9 ) lS the Laplace exponent of this subordinator. 

Knobloch and Kyprianou [56] show that for appropriate values of c,p > 0, 

M t := £ e-^'X^W^ict + logXiit)), t>0, (1.10) 

ieSf 

is a martingale, where for 9 > 0, y(9) = c9 — <p(9) is the Laplace exponent of a 
spectrally negative Levy process, is the g-scale function with respect to y/ and 
J^ e is the set of indices of fragments at time t whose genealogical line of descent 
has the property that, for all 5 < t, its ancestral fragment at time s is larger than 
e~ cs . One may think of the sum as being over a 'thinned' version of the original 
fragmentation process. That is to say, an adjustment of the original process in which 
blocks are removed if, at time t > 0, they become smaller than e~ cr . Removed blocks 
may therefore no longer contribute fragmented mass to the on-going process of 
fragmentation. 

Analysis of this martingale in [56], in particular making use of known properties 
of scale functions, allows the authors to deduce that, under mild conditions, there 
exists a unique constant p* > such that whenever c > <j>'(p*) 

sup-^^V), (1.11) 

as 1 1 00 on the event that the index set remains non-empty for all t > (i.e. the 
thinned process survives), which itself occurs with positive probability. This result is 
of interest when one compares it against the growth of the largest fragment without 
restriction of the index set. In that case it is known that 
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logX ( (0 



-►*V). 



t 



as 1 1 00 • Intuitively speaking (1.11) says that the thinned fragmentation process will 
either become eradicated (all blocks are removed in a finite time) or the original 
fragmentation process will survive the thinning procedure and the decay rate of the 
largest block is unaffected. 

A more elaborate martingale, similar in nature to (1.10) and also built using a 
scale function, was used by Krell [62] for a more detailed analysis of different rates 
of fragmentation decay that can occur in the process. In particular, a fractal analysis 
of the different decay rates is possible. 

Levy processes conditioned to stay positive: What is important in many of the 
examples above is the behaviour of the underlying spectrally negative Levy process 
as it exists a half line. In contrast, one may consider the behaviour of such Levy 
processes conditioned never to exit a half line. This is of equal practical value from 
a modelling point of view, as well as having the additional curiosity that the condi- 
tioning event can occur with probability zero. 

Take, as usual, X to be a spectrally negative Levy process and recall the definition 
(1.6) of Tq. Assume that </(0+) > 0, then it is known that for all t > and x,y > 0, 



where e is an independent and exponentially distributed random variable with unit 
mean, defines the semigroup of a conservative strong Markov process which can be 
meaningfully called a spectrally negative Levy process conditioned to stay positive. 
Moreover, it turns out to be the case that the laws {Pi : x > 0} can be described 
through the laws {P x : x > 0} via a classical Doob /z-transform. Indeed, for all A e & t 
and x > 0, 



Hence a significant amount of probabilistic information concerning this condition- 
ing is captured by the scale function. See Chaumont and Doney [27] for a complete 
overview. In a similar spirit Chaumont [25, 26] also shows that scale functions can 
be used to describe the law of a Levy process conditioned to hit the origin con- 
tinuously in a finite time. Later on in this review we shall see another example of 
conditioning spectrally negative Levy processes due to Lambert [73] which again 
involves the scale function in a similar spirit. 



P}{X t e dy) = \imP x (X t € dy, t < e/q\T > e/q), 



(1.12) 




Chapter 2 

The general theory of scale functions 



2.1 Some additional facts about spectrally negative Levy 
processes 

Let us return briefly to the family of spectrally negative Levy processes and remind 
ourselves of some additional deeper properties thereof which, as we shall see, are 
closely intertwined with the properties of scale functions. 

Path variation. Over and above the assumption of spectral negativity in (1.1), the 
conditions 



are necessary and sufficient for X to have paths of bounded variation. In that case it 
may necessarily be decomposed uniquely in the form 



where 8 > and {S t : t > 0} is a pure jump subordinator. 

Regularity of the half line: Recall that irregularity of (— °°,0) for for X means 
that P(t q > 0) = 1 where 



Thanks to Blumenthal's zero-one law, the only alternative to irregularity of (— °°,0) 
for is regularity of (— °°, 0) for meaning that P(Tq > 0) = 0. Roughly speaking 
one may say in words that a spectrally negative Levy process enters the open lower 
half line a.s. immediately from if and only if it has paths of unbounded variation, 
otherwise it takes an a.s. positive amount of time before visiting the open lower half 
line. In contrast, for all spectrally negative Levy processes we have P(t^ = 0) = 1, 
where 




X, = 8t-S„t >0, 



(2.1) 



T " = inf{f > : X, < 0}. 



T+ = inf{f > : X, > 0}. 
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That is to say, there is always regularity of (0,°°) for 0. 

Creeping. When a spectrally negative Levy process issued from x > enters 
(—oo,0) for the first time it may do so either by a jump or continuously. That is 
to say, on the event jTrT < oo) either X,- — or X_- < 0. The former behaviour 
is called creeping downwards (with positive probability). We are deliberately vague 
here about the initial value x > under which creeping occurs as it turns out that if 
Py(X x - = 0) > for some y > 0, then P X (X T - = 0) > for all x> 0. It is known that 
the only spectrally negative Levy processes which can creep downwards are those 
processes which have a Gaussian component. Note that, thanks to the fact that there 
are no positive jumps, a spectrally negative Levy process will always a.s. creep up- 
wards over any level above its initial position, providing the path reaches that level. 
That is to say, for all a > x, P X (X T + — a; T+ < oo) = P x (t+ < oo) where 

T + = inf{f > : X, > a}. 

Exponential change of measure. The Laplace exponent also provides a natural 
instrument with which one may construct an 'exponential change of measure' on 
X in the sprit of the classical Cameron-Martin-Girsanov transformation, a tech- 
nique which plays a central role throughout this text. The equality (1.2) allows for a 
Girsanov-type change of measure to be defined, namely via 

= Q c{X,-x)- W (c)t ^ t > q (2 2) 

for any c such that | y/(c) | < oo. Note that that it is a straightforward exercise to show 
that the right hand side above is a martingale thanks to the fact that X has station- 
ary independent increments together with (1.2). Moreover, the absolute continuity 
implies that under this change of measure, X remains within the class of spectrally 
negative processes and the Laplace exponent of X under /* is given by 

\lf c (6) = y(d+c)-\if{c), (2.3) 

for > — c. If n c denotes the Levy measure of (X,P C ), then it is straightforward to 
deduce from (2.3) that 

n c (6x) =e cx n(dx), 

for x < 0. 

The Wiener-Hopf factorization. Suppose that X, := mf s < t X s and e q is an inde- 
pendent and exponentially distributed random variable with rate q > 0, then the 
variables X tq — and — are independent. The fundamental but simple concept 
of duality for Levy processes, which follows as a direct consequence of stationary 
and independent increments and cadlag paths states that {X t —X( t _ s y : s < t} is 
equal in law to {X s : s < t }. This in turn implies that X t — X, is equal in distribution 
to X, := sup s<t X s . It follows that for all fldi, 



d^ 
dP x 
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E (e ie ^ ) = E (e ie ^ )E (e ie - e <? ) . (2.4) 

This identity is known as the Wiener-Hopf factorization. Note that the reasoning 
thus far applies to the case that X is a general Levy process. However, when X is a 
spectrally negative Levy process, it is possible to express the right hand side above 
in a more convenient analytical form. Let a > 0, since t A T+ is a bounded stopping 
time and X tAz + < a, it follows from Doob's Optional Stopping Theorem applied to 
the exponential martingale in (2.2) that 

£ ( e *<«>W-«(' A tf>)=l. 

By dominated convergence and the fact that X z + = a on T+ < °° we have, 

P(X e? >a)=P(T+<e 9 ) 

= £(e-^ + V<»)) 

= e -*^ fl . (2.5) 

The conclusion of the above series of equalities is that the quantity X tq is exponen- 
tially distributed with parameter <P(q) and therefore for all 0ef, 

E(^) = -^-. (2.6) 

It follows from the factorization (2.4) that for all 9 e M 

£(e *^) = -g-*(g>- ifl . (2.7) 

As we shall see later on, a number of features described above concerning the 
Wiener-Hopf factorization, including the previous identity, play an instrumental role 
in the theory and application of scale functions. 



2.2 Existence of scale functions 

We are now ready to show the existence of scale functions. That is to say, we will 
show that for Definition 1.4 to make sense, the function (y/(/3) — q)~ l is genuinely 
a Laplace transform in j3 . This will largely be done through the Wiener-Hopf fac- 
torization in combination with an exponential change of measure. 

Theorem 2.1. For all spectrally negative Levy processes, q-scale functions exist for 
all q>0. 

Proof. First assume that y/(0+) > 0. With a pre-emptive choice of notation, define 
the function 
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W W = ^+) P ^>°)- (2-8) 

Clearly W(x) —0 for x < 0, it is right-continuous since it is also proportional to the 
distribution function P(— < x). Integration by parts shows that, on the one hand, 

= (2 ' 9) 

On the other hand, taking limits as q 1 in (2.7) gives us the identity 

E(e WK ~) = -y/(0+) If " 



where, thanks to the assumption that y/(0+) > 0, we have used the fact that <J>(0) = 
and hence 

lim « = lim = lim = y/(0+). 

q\0 <P(q) glO <P(q) 0J.O 9 

A straightforward argument using analytical extension now gives us the identity 

^)-^- 

Combining (2.9) with (2.10) gives us (1.4) as required for the case q = and 
i/(0+) >0. 

Next we deal with the case where q > or q = and y/(0+) < 0. To this end, 
again making use of a pre-emptive choice of notation, let us define 

ff(«>(i)=e^% w (i), (2.11) 

where plays the role of W but for the process (X,P ^). Note in particular 

that by (2.3) the latter process has Laplace exponent 

v*( 4 )(0) = v(e +*(«)) -9. (2A ^ 

for > 0, and hence V*( ? ) (0+) = W'i'&il)) > 0, which ensures that is well 

defined by the previous part of the proof. Taking Laplace transforms we have for 
P><P(q), 
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j~ e-0V«)(jc)djc - ^ ^- 0{q))x W 0(q) {x)Ax 

1 



1 



thus completing the proof for the case g > or g = and y/(0+) < 0. 

Finally, the case that q = and y/(0+) = can be dealt with as follows. Since 
W(p( q ){x) is an increasing function, we may also treat it as a distribution function of 
a measure which we also, as an abuse of notation, call W^ q y Integrating by parts 
thus gives us for j3 > 0, 



[o. 



1 "'" B, *' (d ' ) = ^)- (2 - 13) 



Note that the assumption </(0+) = 0, implies that <P(0) = 0, and hence for 9 > 0, 
limy*( q) (9) =lim[xi/(9 + 0(q))-q] = y/(6). 

One may appeal to the Extended Continuity Theorem for Laplace Transforms, see 
Feller (1971) Theorem XIII. 1.2a, and (2.13) to deduce that since 

km/ e-^W^)(dx) = -^-, 

then there exists a measure W* such that W* (x) := W* [0,x] = lim^o W,p^ (x) and 

Integration by parts shows that the right-continuous distribution, 

W(x) :=W*[0,x] 



satisfies 



for j3 > as required. □ 

Let us return to Theorem 1 .2 and show that, now that we have a slightly clearer 
understanding of what a scale function is, a straightforward proof of the aforemen- 
tioned identity can be given. 

Proof (of Theorem 1.2). First we deal with the case that q = and y/(0+) > as in 
the previous proof. Since we have identified W(x) = PxQL^ > 0)/y/(0+), a simple 
argument using the law of total probability and the Strong Markov Property now 
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yields forx e [0,a] 

Px(£~ > 0) 

= E x (p x (X~>0\&r+)) 

= E.^^PM- > 0)) +E x (i (t , >Tq) P Xt _ (£. > 0)) . (2.14) 
The first term on the right hand side above is equal to 

P o (X»>0)P,(t+ <T "). 

The second term on the right hand side of (2.14) is more complicated to handle 
but none the less is always equal to zero. To see why, first suppose that X has no 
Gaussian component. In that case it cannot creep downwards and hence X T - < 

and the claim follows by virtue of the fact that P X (X X > 0) = for x < 0. If, on the 
other hand, X has a Gaussian component, the previous argument still applies on the 

event {X,- < 01, however we must also consider the event X,- = 0. In that case we 

o T o 

note that regularity of (— °°,0) for foiX implies that PiX^ > 0) = 0. Returning to 
(2.14) we may now deduce that 

and clearly the same equality holds even when x < as both left and right hand side 
are identically equal to zero. 

Next we deal with the case q > 0. Making use of (1.2) in a, by now, familiar way, 
and recalling that X T + — a, we have that 

= e-*W(—)ff W (T+<T -) 

= c -*(Q)(a-x) W *iQ)( X ) 

ffM(i) 

Finally, to deal with the case that q = and </(0+) < 0, one needs only to take 
limits as q I in the above identity, making use of monotone convergence on the 
left hand side and continuity in q on the right hand side thanks to the Continuity 
Theorem for Laplace transforms. □ 

Before moving on to looking at the intimate relation between scale functions 
and the so-called excursion measure associated with X, let us make some further 
remarks about the scale function and its relation to the Wiener-Hopf factorization. 
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It is clear from the proof of Theorem 2.1 that the very definition of a scale func- 
tion is embedded in the Wiener-Hopf factorization. The following corollary to the 
aforementioned theorem reinforces this idea. 

Corollary 2.2. For x>0, 

P(-X eq edx) = -^-^W {q \dx)-qW {q \x)dx. (2.16) 

Proof. A straightforward argument using analytical extension allows us to re-write 
the Wiener-Hopf factor (2.7) as a Laplace transform, 

for j3 > 0. Note in particular, when j3 = <P(q) the right hand side is understood in 
the limiting sense using LHopital's rule. The identity (2.16) now follows directly 
from the Laplace transform (1.4) and an integration by parts. □ 



2.3 Scale functions and the excursion measure 

Section 2.2 shows that there is an intimate relationship between scale functions and 
the Wiener-Hopf factorization. The Wiener-Hopf factorization can itself be seen as a 
distributional consequence of a decomposition of the path of any Levy process into 
excursions from its running maximum, so-called excursion theory; see for example 
the presentation in Greenwood and Pitman [45], Bertoin [14] and Kyprianou [67]. 
Let us briefly spend some time reviewing the theory of excursions within the context 
of spectrally negative Levy processes and thereafter we will show the connection 
between scale functions and a key object which plays a central role in the latter 
known as the excursion measure. 

The basic idea of excursion theory for a Levy process is to provide a structured 
mathematical description of the successive sections of its trajectory which make 
up sojourns, or excursions, from its previous maximum. In order to do this, we 
need to introduce a new time-scale which will help us locate the original times at 
which X creates new maxima. Specifically we are interested in a (random) function 
of original time, or clock, say L = {L t : t > 0}, such that L increases precisely at 
times when X creates a new maximum; namely the times {t : X t = X t }. Moreover, 
L should have a regenerative property in that if T is any F-stopping time such that 
Xt = Xj on {T < °°} then {(XT+t —Xt,Lt+i — Lt) : t > 0} has the same law as 
{(X t ,L t ) : t > 0} under P. The process L is referred to as local time at the maximum. 

It turns out that for spectrally negative Levy processes, there is a natural choice 
of L which is simply L — X. Indeed it is straightforward to verity that X has the 
required regenerative property. We shall henceforth proceed with this choice of L. 
Now define 

L~ x = inf{s > : L s > t}. 
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Note that in the above definition, t is a local time and is the real time which is 
required to accumulate t units of local time. On account of the fact that L = X, it is 
also the case that Ly l is the first passage time of the process X above level t. In the 
notation of the previous section, LJ X = T, + . 

Standard theory tells us that, for all t > 0, both LJ l and Lj} are stopping times. 
Moreover, it is quite clear that whenever ALJ 1 := L~} — LJ X > 0, the process X 
experiences an excursion from from its previous maximum X. For such local times 
t > the associated excursion shall be written as 

E t = {£,{s) :=X ir i -X L - l+s :0<s< AL; 1 }. 

For local times t > such that ALJ X — we define e t — d where d is some isolated 
state. For each t > 0, £ t takes values in the space of excursions, S ', the space of real 
valued right continuous left limited paths killed at the first hitting time of (— °°,0]. 
We suppose that § is endowed with the sigma-algebra generated by the coordinate 
maps, say ^. 

The key feature of excursion theory in the current setting is that {(t,£ t ) : t > 
and e, ^ d} is a Poisson point process with values in (0,°o) x £ , with intensity 
measure df x dn, where n is a measure on the space (S 

Let us now return to the relationship between excursion theory and the scale 
function. First write £ for the lifetime and e for the height of the generic excursion 
e e S. That is to say, for e e 

£ = mi{s > : e(s) < 0} and e = sup{e(s) : s < Q. 

Choose a > x > . Recall that L = X we have that L + —X+ = a —x. Using this 

a— x a—x 

it is not difficult to see that 

{X x + > —x} = {Vf < a -x and e, ^ d, £, < t +x}. 

Let be the Poisson random measure associated with the Poisson point process of 
excursions; that is to say for all dt x n(de) -measurable sets B e [0,°°) x S, N(B) = 
card{(f,£ ( ) £ B,t > 0}. With the help of (2.15) we have that 

>o) 

P(Vf <a-x and e, ^ d,e, < t+x) 
P(N(A)=0), 

where A — {(t , S t ) : t < a — x and e, > t +x}. Since N(A) is Poisson distributed with 
parameter / 1 A dt n(de) = f£~ x n(e > t + x) At = J"n(e > t)dt we have that 



W(x) 
W(a) 




(2.17) 



W(x) 
W{a) 
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This is a fundamental representation of the scale function W which we shall re- 
turn to in later sections. However, for now, it leads us immediately to the following 
analytical conclusion. 

Lemma 2.3. For any q > 0, the scale function is continuous, almost every- 
where differentiate and strictly increasing. 

Proof. Assume that q = 0. Since a may be chosen arbitrarily large, continuity and 
strict monotonicity follow from (2.17). Moreover, we also see from (2.17) that the 
left and right first derivatives exist and are given by 

W'_{x) =n(e >x)W(x) and W[(x) = n(e > x)W(x). (2.18) 

Since there can be at most a countable number of points for which the monotone 
function n(e > x) has discontinuities, it follows that W is almost everywhere differ- 
entiate. 

From the relation (2.1 1) we know that 

ffM(x) = e*«% (?) W, (2.19) 

and hence the properties of continuity, almost everywhere differentiability and strict 
monotonicity carry over to the case q > 0. □ 

Let us write as an abuse of notation ffW e C 1 (0, °°) to mean that the restriction of 
to (0,°o) belongs to C 1 (0,°°). Then the proof of the previous lemma indicates 
that W G C (0,°°) as soon as the measure n(e € dx) has no atoms. It is possible 
to combine this observation together with other analytical and probabilistic facts to 
recover necessary and sufficient conditions for first degree smoothness of W. 

Lemma 2.4. For each q>0, the scale function ff™ belongs to C 1 (0, °°) if and only 
if at least one of the following three criteria holds, 

(H) /(-i,o) \x\n(dx)=<*> 
(Hi) TI(x) :— °°, —x) is continuous. 

Proof. Firstly note that it suffices to consider the case that q — and 1/(0+) > 0, 
i.e. the case that the Levy process oscillates or drifts to +°°. Indeed, in the case that 
y/(0+) < or q > 0, one recalls that <P(0) > 0, respectively <P(q) > 0, and we 
appeal to (2.19), (2.11) where W<p(o)> respectively W$( ? ), is the scale function of a 
spectrally negative Levy process which drifts to +°o. 

If it is the case that there exists an x > 0, such that n(e — x) > 0, then one may 
consider the probability law n(-\s = x) = «(•,£ = x)/n(e = x). Assume now that 
X has paths of unbounded variation, that is to say (i) or (ii) holds, and recall that 
in this case is regular for (— °°,0). If we apply the strong Markov property under 
n(-\e=x) to the excursion at the time T x :=inf{f > : e(t) =x), the aforementioned 
regularity of X implies that (jc,°°) is regular for {x} for the process e and therefore 
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that e > x, under n(-\e=x). However this constitutes a contradiction. It thus follows 
that n(e = x) = for all x > 0. 

Now assume that X has paths of bounded variation. It is known from [83] and 
Proposition 5 in [106] that when X is of bounded variation the excursion measure of 
X reflected at its supremum begins with a jump and can be described by the formula 

n(F(e(s),0 <s < 0) = i J° TI{dx)E- x (F{X s ,0 < s < T Q )) , 

where F is any nonnegative measurable functional on the space of cadlag paths, E x 
denotes the law of the dual Levy process X = —X and T~ = inf{s > : X s < x}, 
i£i Recall that T+ = inf j> > : X s > z}, z € M. Hence, it follows that 

n{e>z) = \ [ n(dx)P- x [ sup X s > z | 

8 ■/(— .0) \0<s<z- J 

= in ( -oo,- z)+ i^ o n(cbc)p_,(T+<T -) 
= ln(-^-z) + y [zQ n( dx )p^ +x <^) 

= ^(-00,-2) + ^ ^ H(dx)P (ti,_ z < T+ ) 

-s n (--"4jL n « fc) ( 1 - ! w)' 

where in the last equality we have used Theorem 1 .2 with q = 0. From this it follows 
that 

n(e = z) = - 5 n({-z})^J ) , 

and hence the required property for that n(e = z) =0 (which in turn leads to the 
C l (0,°°) property of the scale function) follows as soon as II has no atoms. That is 
to say n is continuous. □ 



The proof of the above lemma in the bounded variation case gives us a little more 
information about non-differentiability in the case that II has atoms. 

Corollary 2.5. When X has paths of bounded variation the scale function does 
not possess a derivative at x > (for all q>0) if and only if II has an atom at —x. 
In particular, ifTi has a finite number of atoms supported by the set {—x\ , • • • , — x n } 
then, for all q> 0, W (<?) G C 1 ((0,°°)\{jci, • • • ,*„}). 

We conclude by noting that the above results give us our first analytical impres- 
sions on the 'shape' and smoothness of In Section 3 we shall explore these 
properties in greater detail. 
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2.4 Scale functions and the descending ladder height process 

In a similar spirit to the previous section, it is also possible to construct an excursion 
theory for spectrally negative Levy processes away from their infimum. Indeed it 
is well known, and easy to prove, that the process X reflected in its past infimum 
X — X := {X, —X_[ : t > 0}, where X^ := inf^X,, is a strong Markov process with 
state space [0,°°). Following standard theory of Markov local times (cf. Chapter IV 
of [14]), it is possible to construct a local time at zero forX — X_ which we henceforth 
refer to as L = {L t : t > 0}. Its right-continuous inverse process, Lr x : = {L^ 1 : t > 0} 
where LJ X — inf{s > : L s > f}, is a (possibly killed) subordinator. Sampling X 
at L _1 we recover the points of local minima of X. If we define H, = X-£-\ when 

L^ 1 < oo, with H t = oo otherwise, then it is known that the process H = {H t : t > 0} 
is a (possibly killed) subordinator. The latter is known as the descending ladder 
height process. 

The starting point for the relationship between the descending ladder height pro- 
cess and scale functions is given by the following factorization of the Laplace expo- 
nent of X, 

iKe) = (e-0(o))<H0), e>o, (2.20) 

where <j> , is the Laplace exponent of a (possibly killed) subordinator. See e.g. Section 
6.5.2 in [67]. This can be seen by recalling that <£(()) is a root of the equation 
y(0) = and then factoring out (9 — 0(0)) from y/ by making some integration by 
parts to deduce that the term <j>(6) must necessarily take the form 

<j>{d) = K + ^d+ f (l-e- ex )Y(dx), (2.21) 

y(o,oo) 

where K.E, > and Y is a measure concentrated on (0,°°) which satisfies f^ ^(1 A 
x)Y(dx) < oo. Indeed, it turns out that 

e -*(°)»n(-oo,-«)da forx>0, (2.22) 

t, =o 2 /2 and K = E(X\ ) V = i//(0+) V 0. 

Ultimately the factorization (2.20) can also be derived by a procedure of analyti- 
cal extension and taking limits as q tends to zero in (2.4), simultaneously making use 
of the identities (2.6) and (2.7). A deeper look at this derivation yields the additional 
information that is the Laplace exponent of H; namely 0(0) = — \ogE(e~ eH '). 

In the special case that <I>(0) = 0, that is to say, the process X does not drift to -oo 
or equivalently that y/(0+) > 0, it can be shown that the scale function W describes 
the potential measure of H. Indeed, recall that the potential measure of H is defined 
by 

/ dt-P(H,edx), forx>0. (2.23) 
Jo 

Calculating its Laplace transform we get the identity 
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f ( At ■ P(H t e dx)e 



1 



for 9 > 0. 



Jo Jo 



(2.24) 



Inverting the Laplace transform on the left hand side we get the identity 



W(x) = [ &t-P(H, € [0,x]), 



x>0. 



(2.25) 



It can be shown similarly that in general, when 4>(0) > 0, the scale function is 
related to the potential measure of H by the formula 



The connection between W and the potential measure of H turns out to be of 
importance at several points later on in this exposition. 



2.5 A suite of fluctuation identities 

Let us now expose the robustness of scale functions as a natural family of functions 
to describe a whole suite of fluctuation identities concerning first and last passage 
problems. We do not offer full proofs, concentrating instead on the basic ideas that 
drive the results. These are largely concentrated around the Strong Markov Property 
and earlier established results on the law of X e and — X~ . 

Many of the results in this section are taken from the recent work [15, 84, 9, 83]. 
However some of the identities can also be found in the compound Poisson setting 
from earlier Soviet-Ukranian work; see for example [57, 58, 59, 22, 61, 60, 99] to 
name but a few. 



2.5.1 First passage problems 

Recall that for all a e K, T+ = inf{f > : X t > a} and T ~ = inf{f > : X t < 0}. We 
also introduce for each q > 0, 




x>0. 



(2.26) 




(2.27) 



Theorem 2.6 (One- and two-sided exit formulae). 



C i) For any i€i and q>0 




(2.28) 
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where we understand q / (q) in the limiting sense for q — 7 so that 

p(r -^ , Jl-y/(0+)W(*)ify/(0+)>0 

P ^o = ifv /(0 + )<0- (2 ' 29) 

( ii) For allx>0 and q>0 

2 

E x (e-<*°l {x _ =0] ) = ^{ W ^'(x)-0(q)W^(x)}, (2.30) 

where the right hand side is understood to be identically zero if a — and oth- 
erwise the derivative is well defined thanks to Lemma 2.4. 
( Hi) For any x<a and q>0, 

^(e--o 1(T _ <Ta+) )=^) W -^)(a)^M. (2.31) 



Sketch proof, (i) Making use of Corollary 2.2 we have for q > and x > 0, 

^( e " 9T °" 1 (T -<oo ) ) =Px{e q >To) 

= Px(* eq < 0) 
= P{-X* q >x) 
= l-P(-X eq <x) 

= l+ q J*W(<!)(y)d y -^-WM(x) 

= Z^(x)-^ n W M (x). (2.32) 
sP (q) 

The proof for the case q = follows by taking limits as q I on both sides of the 
final equality in (2.32). 

(iii) Fix q > 0. We have for x > 0, 

^(e^l^ <tf p=^(e-^l (ljr<-) )-^(e-^ri (tf<1ir) ). 

Applying the Strong Markov Property at T+ and using the fact that X creeps up- 
wards, we also have that 

Piecing the previous two equalities together and appealing to (2.28) and (1.5) yields 
the desired conclusion. The case that q — is again handled by taking limits as q 1 
on both sides of (2.31). Here we have used the discussion in Section 2.4. 

(ii) Suppose that q = and </(0+) > (equivalently the descending ladder 
height process is not killed) then the claimed identity reads 
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P x (X T -=0) = —W'(x), x>0. 



(2.33) 



Note that the probability on the left hand side is also equal to the probability that 
the descending ladder height subordinator creeps over x, P(Hf = x), where T x = 

inf{f >0:H t > x}. The identity (2.33) now follows directly from a classic result 
of Kesten [53] which shows that the probability that a subordinator creeps is non- 
zero if and only its drift coefficient is strictly positive, in which case it is equal to 
the drift coefficient multiplied by its potential density, which necessarily exists. It 
follows together, with (2.22) and (2.26), that P(H fx =x) = o 2 W'{x) /2. 

When q > or y/(0+) < 0, the formula is proved using the change of measure 
(2.2) with c = <l>(q), q>0, and the above result together with (2.1 1). □ 

Note that part (ii) above tallies with the earlier mentioned fact that spectrally neg- 
ative Levy processes do not creep downwards unless they have a Gaussian compo- 
nent. Note also that when ff^O Lemma 2.4 tells us that the derivative appearing on 
the right hand side of the density is everywhere defined for x > 0. 

We also give expressions for the expected occupation measure of X in a given 
Borel set over its entire lifetime as well as when time is restricted up to the first 
passage times T+, Tq and 



Such expected occupation measures are generally referred to as resolvents and play 
an important role in establishing, for example, deeper identities concerning first 
passage problems such as the one presented in Theorem 1.3. 

Theorem 2.7 (Resolvents). 

C i) For all a > x > 0, q > and Borel set AC [0, a], 



(wM(x)W(«\a-y) 



Jo 6 w^j-y^- 



( ii) For all a>x and Borel set A C (— °o,a], 



■W ( - q \x-y)\dy. 



{x,eA,t<tt} 



df 



: j | e -*(<7)( a -*) (a-y)-W W (x - y) } dy. 



( Hi) For allx>0 and Borel set A C [0, °°), 



j\-n {XteAj<T - } dt] = j A {c-^yW^(x)-W^(x-y)}dy. 



(iv) For all x gl and Borel set A C. 



jfV*V €A} d/| = J a {&(q)e-*byy- W<-*\-y)}dy. 
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Sketch proof. We give an outline of the proof of (iii) from which the proof of (i) 
easily follows. The remaining two identities can be obtained by taking limits of the 
barriers in (i) relative to the initial position. As usual we shall perform the relevant 
analysis in the case that q > 0. The case that q = follows by taking limits as q 1 0. 
We start by noting that for all x,y > and q > 0, 

R^(x,dy) := I e-* P x (X t e dy,T > t)dt = -P x (X eq e dy,X e? > 0), 

where is an independent, exponentially distributed random variable with param- 
eter q > 0. 

Appealing to the Wiener-Hopf factorization, specifically that X eq — X eq is inde- 
pendent of X^ q , and that X eq — X eq is equal in distribution to X eq , we have that 

*<«>(*, dy) = -P((X eq -X eq )+X eq € dy-x,-X Bq <x) 

= \j m P ( G ^) P (*e ? G dy - x + z ) 1 . 

Recall however, that X tq is exponentially distributed with parameter <P(q). In addi- 
tion, the law of ~X eq has been identified in Corollary 2.2. Putting the pieces together 
and making some elementary manipulations the identity in (iii) follows. 

Now suppose we denote the left hand side of the identity in (i) by U^ q \x,A). 
With the help of the Strong Markov Property we have that 

qU^(x,dy) = P x (X eq € dy,Xe q > 0,X eq < a) 
= P x (X eq edy,X eq >0) 

-P x (X eq e dy,X eq > 0,X eq > a) 
= P x (X eq edy,X eq >0) 

-P X (X X =a,T< e q )P a (X eq e dy,X tq > 0). 

The first and third of the three probabilities on the right-hand side above have been 
computed in the previous paragraph, the second probability may be written 

The result now follows by assembling the relevant pieces. □ 

On a final note, the above resolvents easily lead to further identities of the type 
given in Theorem 1.3. Indeed, suppose that N is the Poisson random measure as- 
sociated with the jumps of X. That is to say, N is a Poisson random measure on 
[0,°°) x (-°o,0) with intensity df x IT(dx). Recall that T := T+ A Xq . Then with the 
help of the Compensation Formula, we have that for x G [0,a], A any Borel set in 
[0,a) and B any Borel set in (— °o,0) and q > 0, 
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P x (e- qT ;X r GB,X r - eA) 

= Ex (/oh /— ,0) e_, ' 1 ^-<«&->0A-eA) 1 (y€i>-*-)^(d* x d ?)) 
= E^V% <T) I7(B-X f )l(x, e A)d/) 

= / II(B -?)£/(«> (jc,dy), (2.34) 



7A 

where, as noted earlier, 



/■CO 

£/(«) (jc, dy) = / e-^(X, e dy, T > f )dr . 
Jo 

Observe that the fact that B is a Borel subset of (— °°, 0) , allow us not to consider the 
event where the process leaves the interval from below by creeping. Nevertheless, 
the probability of this event has been calculated in (2.30). 



2.5.2 First passage problems for reflected processes 

The list of fluctuation identities continues on when one considers the reflected pro- 
cesses F := {X, -X, : t > 0} and Y := {X, -X, : t > 0}. Note that it is easy to prove 
that both of these processes are non-negative strong Markov processes. We shall 
henceforth denote their probabilities by {P x :x E [0, °°)} and {P x :i£ [0, °°)}. Note 
that (Y,P X ) is equal in law to {(xVX t — X t ) : t > 0} under P and (Y_,P X ) is equal in 
law to {X t — (X_ t A —x) : t > 0} under P. The following theorem is a compilation of 
results taken from [9] and [83]. We do not offer proofs but instead we shall settle for 
remarking that the proofs use similar techniques to those of the previous section. 
For convenience, let us define for a > 0, 

g_ a = inf{f > : T, > a} and a a — inf{f > : Y, > a}, 

where Y_ t =X t -X t and Y, = X, -X,. 

Theorem 2.8. Suppose that a > 0, x £ [0, a] and q>0. Then 

(i) E x (e-^)=Z^(x)/Z^(a), 

(ii) taking W^' (a) is the right derivative ofW^ at a, 

E x {e-«°«)=Z l «\a-x)-qW ( «> (a-x)W W (a) /W^' (a), 
(Hi) for any Borel set A e [0,a), 

j-ao 



zWa) 



W {q) (a-y)-W {q) (x~y) L dy, 
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and 

(iv) for any Borel set A € [0,a), 



Jo 



1 {Y,eA,t<a a } dt 



Ja\ Wi q) '(a) 
+ f lwl<Hx-a)™^ 



Chapter 3 

Further analytical properties of scale functions 



3.1 Behaviour at and +°° 



Ultimately we are interested in describing the 'shape' of scale functions. We start 
by looking at their behaviour at the origin and +°°. In order to state the results 
more precisely, we recall from (2.1) that when X has paths of bounded variation, 
we may write it in the form X t = St — S t , t > 0, where 8 > and S is a pure jump 
subordinator. 

Lemma 3.1. For all q>0, W^ q \0) — if and only ifX has unbounded variation. 
Otherwise, whenX has bounded variation, W^ q \0) = 1/8. 

Proof. Note that for all q > 0, 

W ( «'(0) = lim Hp e-P x W M (x)dx 
/3W0 

B 

= Hm 

= lim-^. (3.1) 

Recall the spatial Wiener-Hopf factorization of y/ in (2.20), 

W(d) = (d-0(O))<l>(d), e>o, 

and that denotes the Laplace exponent of the downward ladder height subordinator 
H. It follows that 

wW(0)=lim-^-= lim 

Now, observe that limg^oo (j8) < °°, if and only if the Levy measure of H is a 
finite measure and its drift is 0. We know this happens if and only if X has paths 
of bounded variation, as we already mention that this is the only case in which 
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is irregular for (— oo,0) and hence starting from it takes a strictly positive amount 
of time to enter the open lower half line. The first claim of the Theorem follows. 
Now, assume thatX has paths of bounded variation. In this case, one may use (2.1) 
to write more simply 



wiP) = sp-f (l-e-^)n(dx). 

J(0,oo) 



An integration by parts tells us that 



r-oo 

/ e-"*n(-oo,-jc)dx, 
Jo 



(3.2) 



and hence, noting that bounded variation paths (equivalently j |jc| JT(dbc) < °°) 
necessarily implies that ^ JT(— oo, — x )dx < oo, it follows that yf(j3)/j3 -> 8 as 

j3 f °°- From (3.1) we now see that, for all q > 0, ffW (0) = 1/5 as claimed. 

Alternatively, one can prove the result making an integration by parts for y/ in- 
stead of using the spatial Wiener-Hopf factorization. □ 

Returning to (2.15) we see that the conclusion of the previous lemma indicates 
that, precisely when X has bounded variation, 

^o(T a + <T -)^^||>0. (3.3) 

Note that the stopping time Tq is defined with strict entry into (— °°,0). Hence when 
X has the property that is irregular for (— °°,0), it takes an almost surely positive 
amount of time to exit the half line [0,°°). Since the aforementioned irregularity is 
equivalent to bounded variation for this class of Levy processes, we see that (3.3) 
makes sense. 

Next we turn to the behaviour of the gradient of at 0. 



Lemma 3.2. For all q>0 we have 
Wr(«)'( +) 



2/ a 2 , when a ^ or J7(-oo,0) = °° 

(n(-°o,0)+g)/(5 2 when a = and J7(-oo,0) < oo, 



where we understand the first case to be +°° when a = 0. 
Proof. Integrating (1.4) by parts we find that for 9 > <P(q), 

W W (0) + f e- 0x W W ' (jc)dx = — ^ . (3 .4) 

Jo \ir{v)-q 

Now assume that X has paths of unbounded variation so that W^(0) = 0. Next 
noting from Lemma 2.3 that a right derivative of at zero always exists, we have 
for each q > 0, 
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e 2 



'(0+) = lim f°° 9e- ex W M '(x)dx = lim 



y(9)- q - 



Then dividing (2.20) by 2 it is easy to proved using (2.22) that the limit above 
is equal to 2/a 2 as required. The expression for W^'(0+) in the first case now 
follows. 

When X has bounded variation, a little more care is needed. Recall that (0+) = 
1/5, we have, 



fff«''(0+) 



= lim 



= lim 



p2 fi -W^(0+)8+W^(0+)foe-^ x n(-oo,-x)6xfj +9j3W^(0+) 

<5j3 - £ Pe-P*n(-°°, -x)dx + q 
1 / °°j3e-^n(-°o,-x)dx + (? 



pt-8 5 J °° e-P*n(-oo, -jc)djc 

n(-°o,o)+ g 

In particular, if JT(— °°,0) = °°, then the right hand side above is equal to °° and 
otherwise, if JT(-oo,0) < °°, then W (<?) '(0+) is finite and equal to (JI(-oo,0) + 
q)/8 2 . □ 

Next we look at the asymptotic behaviour of the scale function at +°o. 

Lemma 3.3. For q>0we have, 

lim e- ^ x W^(x) = l/xi/(0(q)), 

X— S-oo 

and 

limZ^(x)/W ( - q \x)=q/<P(q), 

where the right hand side above is understood in the limiting sense lim^o q/<P(q) = 
0V(l/\i/(0))whenq = 0. 

Proof. For the first part, recall the identity (2.19) which is valid for all q > 0. It 
follows from (2.8) that 

W («) ( X ) = e *M* . 1 pfW^ > 0). (3.5) 

Appealing to (2.12) we note that V / $( (y )(0+) = </(<£(<?)) > (which in particular 
implies that X under P*(«) drifts to +°o) and hence 

lime-*^W^(x) = — 1— -limp* ®(X„ > 0) = * .. . 
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Note that from this proof we also see that W#( ? ) (+°°) = 1/ 

For the second part, it suffices to compare the identity (2.31) as a f °° against 
(2.28). □ 



3.2 Concave-convex properties 

Lemma 3.3 implies that when q > 0, grows in a way that is asymptotically 
exponential and this opens the question as to whether there are any convexity prop- 
erties associated with such scale functions for large values. Numerous applications 
have shown the need to specify more detail about the shape, and ultimately, the 
smoothness of scale functions. See for example the discussion in Chapter 1. In this 
respect, there has been a recent string of articles, each one improving of the last, 
which has investigated concavity-convexity properties of scale functions and which 
are based on the following fundamental observation of Renming Song. Suppose that 
</(0+) >0 and that 

n(x):=n(-co,-x), 

the density of the Levy measure of the descending ladder height process, see (2.22), 
is completely monotone. Amongst other things, this implies that the descending 
ladder height process H belongs to the class of so-called complete subordinators. A 
convenience of this class of subordinators is that the potential measure associated to 
H, which in this case is W, has a derivative which is is completely monotone; see 
for example Song and Vondracek [96]. In particular this implies that W is concave 
(as well as being infinitely smooth). Loeffen [75] pushes this idea further to the case 
of for q > as follows. (Similar ideas can also be developed from the paper of 
Rogers [88]). 

Theorem 3.4. Suppose that —17 has a density which is completely monotone then 
has a density on (0,°°) which is strictly convex. In particular, this implies the 
existence of a constant a* such that is strictly concave on (0,a*) and strictly 
convex on (a*,°°). 

Proof. Recall from (2.11) that we may always write W^(x) = e*W% (9) (i), 
where W<p( 9 ) plays the role of W under the measure P*(*). Recall also from the 
discussion under (2.12) and (2.3) that (X,P*W) drifts to +°° and that the Levy 
measure of X under P*M is given by n 0{q) (dx) = e ^ x n(dx), x e K. It follows 
that under the assumptions of the Theorem that the function n^/ q \ (— °°, —x), x > 0, 
is completely monotone. Hence the discussion preceding the statement of Theo- 
rem 3.4 tells us that W <i> {q){x) exists and is completely monotone. According to the 
definition by Song and Vondracek [96], this makes of W$( q ) a Bernstein function. 

The general theory of Bernstein functions dictates that there necessarily exists a 
triple (a,b,§), where a,b > and t; is a measure concentrated on (0,°o) satisfying 

.[(0H( lAf ^( df ) <°°' such that 
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W 0(q) (x) = a + bx+ [ (l-e- xt )%{dt). 

J (0,00) 

It is now a straightforward exercise to check with the help of the above identity and 
the Dominated Convergence Theorem that for x > 0, 

W {q) "'{x) = f"'(x) + [ (0(q) 3 ^ x - (0(q) - f ) 3 e (*<«>-'>*)§ (df ) 

+ [ (0( ? ) 3 e^ + (f-0(< 7 )) 3 e-^-*(' ? »)^(df), 
•/(*(?),<») 

where f(x) = (a + bx)^*. Hence W (<?) "'(x) > for all x > 0, showing that W iq) ' 
is strictly convex on (0,°°) as required. □ 

Following additional contributions in [71], the final word on concavity-convexity 
currently stands with the following theorem, taken from [78], which overlaps all of 
the aforementioned results. 

Theorem 3.5. Suppose that TI is log-convex. Then for all q > 0, has a log- 
convex first derivative. 

Note that the existence of a log-convex density of — TI implies that TI is log- 
convex and hence the latter is a weaker condition than the former. This is not an 
obvious statement, but it can be proved using elementary analytical arguments. 



3.3 Analyticity in q 

Let us now look at the behaviour of as a function in q. 

Lemma 3.6. For each x > 0, the function q i-> (x) may be analytically extended 
to q € C. 

Proof. For a fixed choice of q > 0, 



.70 



j 1 

V{P)-q 
1 1 

= y(/3)i-?/y(/3) 



1 L^-n^k' w 



for J3 > <P (q) . The latter inequality implies that < qj 1/a(j8 ) < 1 . Next note that 

£ q k W*^ (x) 
k>0 
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converges for each x > where W* k is the k th convolution of W with itself. This is 
easily deduced once one has the estimates 

W< k+ V(x)<^W(x) k+ \ x>0; (3.7) 

which itself can easily be proved by induction. Indeed note that if (3.7) holds for 
k > 0, then by mono tonicity of W, 



W< k+i Hx)< J*-L— W (y) k W(x-y)dy 



(k-l) 

= ^W{x) k+ \ 

Returning to (3.6) we may now apply Fubini's Theorem (justified by the assumption 
that j3 > <P(q)) and deduce that 

/'CO 1 

Jo V ; kk w(B" k+l 



t>, w(PT 

i-oa 

= Yq k e-P x W< k+1 Hx)dx 
k>0 Jo 

r-ao 

= / e~ l}x Yq k W*V c+ V{x)dx. 
JO t>n 



k>0 

Thanks to continuity of W and we have that 

W( q \x) = Y,q k W< k+ V{x),xeW. (3.8) 

<r>0 

Now noting that £jt>o#*W*(* +1 ) (x) converges for all q <G C we may extend the 
definition of for each fixed x > 0, by the equality given in (3.8). □ 

For each c > 0, denote by wj^ the g-scale function for (X,P C ). The previous 
Lemma allows us to establish the following relationship for wj 9 ^ with different val- 
ues of q and c. 

Lemma 3.7. For any geC and c£R such that y/(c) < °° we have 

W^\x)=e cx W c {q ~ v{c) \x), (3.9) 

for all x > 0. 

Proof. For a given eel, such that y/(c) < °o ; the identity (3.9) holds for q— \j/(c) > 
0, on account of both left and right-hand side being continuous functions with the 
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same Laplace transform. By Lemma 3.6 both left- and right-hand side of (3.9) are 
analytic in q for each fixed x > 0. The Identity Theorem for analytic functions thus 
implies that they are equal for all q e C. □ 

Unfortunately a convenient relation such as (3.9) cannot be given for 
Nonetheless we do have the following obvious corollary. 

Corollary 3.8. For each x > the function q Z^' (x) may be analytically extended 
toqeC 

The above results allow one to push some of the identities in Section 2.5 further 
by applying an exponential change of measure. In principle this allows one to gain 
distribution information about the position of the Levy process at first passage. We 
give one example here but the reader can easily explore other possibilities. 

Consider the first passage identity in (2.28). Suppose that v > 0, then | y(v) | < °°, 
and assume it satisfies u > \jf(v) V 0. Then with the help of the aforementioned 
identity together with the change of measure (2.2) we have for ieR, 



-UXr, +VX 

E x e M 



(t <<*>) 



£;( e -(»-v(v))%i (T _ <oo) ) 



where p = u — y/(v) > 0. We can develop the right hand side further using the rela- 
tionship between scale functions given in Lemma 3.7 as well as by noting that 



<P v (p) = sup{A > 
= sup{A > 
= sup{0 >0 
= <P(u) - v. 



y/ v (X)=p} 

y/(A + v) - y/(v) =t/-y(v)} 



Hence we have that 



— UXr, +VX 



(m-Va(v)) [\- vy W (u) (y)dy- - - y(v) e-"ffW(i) ) . (3.10) 
Jo 



<P(u) -v 

Clearly the restriction that u > y/(v) is an unnecessary constraint for both the left and 
right hand side of the above equality to be finite and it would suffice that m, v > 0. In 
particular for the right hand side, when u = y/(v) it follows that <f>(u) = v and hence 
the ratio (u — \j/(v)) / (<P(u) — v) should be understood in the limiting sense. That is 
to say, 

lim ^ = limVA ^T )) -^(0+) = V/>), 

p^0<P v (p) p^Q <P v (p) 

where we have used the fact that <£ v (0) = as 1/^(0+) = y/(v) > 0. 
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Note that the left hand side of (3.10) is analytic for complex u with a strictly 
positive real part. Moreover, thanks to Lemma 3.6 and the fact that 4>(h) is a Laplace 
exponent (cf. the last equality of (2.5)) it is also clear that the right hand side can be 
analytically extended to allow for complex-valued u with strictly positive real part. 
Hence once again the Identity Theorem allows us to extend the equality (3.10) to 
allow for the case that u > 0. The case that u = can be established by taking limits 
as u 1 on both sides. 

Careful inspection of the above argument shows that one may even relax the 
constraint on v > to simply any v such that | y/(v) | < °° in the case that the exponent 
\jf(v) is finite for negative values of v. 



3.4 Spectral gap 

Bertoin [16] showed an important consequence of the analytic nature of scale func- 
tion ffW in its argument q. 

For a > 0, let T := T+ A Tq, then [16] investigates ergodicity properties of the 
spectrally negative Levy process X on exiting [0,a\. The main object of concern is 
the killed transition kernel 

P{x,t,A)=P x {X,eA;t<T). 

Theorem 3.9. Define 

P = in%>0: W { - q) {a) =0}. 

Then p is finite and positive, and for any q < p and x <G (0, a), W^~ q ^ (x) > 0. Fur- 
thermore, the following assertions hold, 

(i) p is a simple root of the entire function q W^~ q \a), 

(ii) The function W^ p ^ is positive on (0,a) and is p-invariant for P(x,t,-) in the 
sense that 

f P(x,t 1 dy)W i - p) (y)^e- pt W { - p) {x),foranyxe (0,a), 

J[0,a] 

(Hi) the measure W^ p ^ (a — x)dx on [0,a] is p-invariant in the sense that 
( dyW { - p) (a-y)P(y,t,dx)=e- pt W { - p) (a-x)dx, 

J[0,a] 

( iv) There is a constant c> such that, for any x £ (0, a), 

\m\& pt P{x 1 t 1 dy)^cW { - p) {x)W ( - p) {a-y)dy 1 

in the sense of weak convergence. 
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A particular consequence of the part (iv) of the previous theorem is that there 
exists a constant d > such that 

P x {T>t)~dW { - p) {x)e-P t . 

as 1 1 °°. The constant p thus describes the rate of decay of the exit probability. By 
analogy with the theory of diffusions confined to compact domains, — p also plays 
the role of the leading eigen-value, or spectral gap, of the infinitesimal generator of 
X constrained to the interval (0,a) with Dirichlet boundary conditions. From part 
(i) of the above theorem, we see that the associated eigen-function is W^ p \ 

Lambert [73] strengthens this analogy with diffusions and showed further that 
for each xG (0,a), 

is a P x martingale which, when used as a Radon-Nikodim density to change measure, 

t 

induces a new probability measure, say, P£. He shows moreover that this measure 
corresponds to the law of X conditioned to remain in (0, a) in the sense that for all 
t >0, 

^(A)=HmP x (A\T>t+s), Ae^, 

The role of the scale function is no less important in other types of related condi- 
tioning. We have already alluded to its relevance to conditioning spectrally negative 
Levy processes to stay positive in the first chapter. Pistorius [83, 82] also shows that 
a similar agenda to the above can be carried out with regard to reflected spectrally 
negative Levy processes. 



3.5 General smoothness and Doney's Conjecture 

Let a > 0, and recall T = T+ A Tq . It is not difficult to show from the identity (1.5) 
that for all q > 0, 

& -i{t^) W {i){x tM ),t>0 

is a martingale. Indeed, we have from the Strong Markov Property that for all x e R 
and t > 

= e -<?(?AT - AT+) W(q) ( X 'AT Q - AT+ ) 

W(D( a ) ' 

where it should be noted that for the second equality we have used the fact that 
{X x +)/W^ (a) = 1 and (X T - ) /W^ (a) = 0. Note in particular, for the last 
equality, X creeps downwards if and only if it has a Gaussian component in which 
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case W^(X T -) = W(?)(0+) = 0, on the event of creeping, and otherwise on the 

event that X T - < 0, it trivially holds that W (<?) (X T - ) = 0. 

Naively speaking this reiterates an idea seen in the previous section that is 
an eigen-function with respect to the infinitesimal generator of X with eigen-value 
q. That is to say, solves the integro-differential equation 

(r-q)W {g) {x) =0on (0,a), 
where F is the infinitesimal generator of X, which is known to be given by 

rf(x) = Hf'(x) + \o 2 f"{x) + / {f(x+y)-f(x) -f'(x)yl {y> _ l} }n(dy), 

J(—<*>fi) 

for all / in its domain. 

The problem with the above heuristic observation is that may not be in the 
domain of F. In particular, for the last equality to have a classical meaning we need 
at least that / G C 2 (0,a) and f,^ j f(x+y)Tl(dy) < °°. This motivates the question 
of how smooth scale functions are. Given the dependency of concavity-convexity 
properties on the Levy measure discussed in Section 3.2 as well as the statement 
of Lemma 2.4, it would seem sensible to believe that a relationship exists between 
the smoothness of the Levy measure and the smoothness of the scale function W^ q \ 
In addition, one would also expect the inclusion of a Gaussian coefficient to have 
some effect on the smoothness of the scale function. Below we give a string of recent 
results which have attempted to address this matter. For notational convenience we 
shall write eC*(0,°°) to mean that the restriction of to (0,°°) belongs to 
theC^O,^) class. 

Theorem 3.10. Suppose that a 2 > 0, then for all q>0, € C 2 (0,°o). 

Proof. As before, it is enough to consider the case where q = and y/(0+) > 0. In 
other case we use the latter with the help of (2.11). Note that it suffices to prove the 
result for q = thanks to (2.11) and (2.12). We know from Lemma 2.4 that when 
a 2 >0,ffW eC l {0,°°). Recall from (2.18) that W'(x) = n(e > x)W(x) and hence 
if the limits exist, then 

W" { x + )=^ W '^- W '^ 
e|0 e 

lim «(e e [x,x + e))W(x) -n(e>x + e)(W(x + e) - W(x)) 
eiO e 

= -lim " (£€ {x ' X + £)) W(x)+n(e>x)W'(x+). (3.11) 
e|0 e 

To show that the limit on the right hand side of (3.11) exists, define a x — inf{/ > 
: £ ( > x} and & t = o(e s : s<t). With the help of the Strong Markov Property for 
the excursion process we may write 
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n(ee [x,x + e)) = n(a x <°°,e ax <x + e,e <x + e) 

= n(l {ax<a0}E(!x<x+e] n{e <x+e|Sf (Tx )) 

= n(l{o x <°°,e ax <x+e} P -e„ x (i:o < T -(x+e))) 

W(e) 



= n(a x < oo, = x) 



W(x + e) 



+« ^<-^<m} w(x + £) J • (3-12) 

We know that spectrally negative Levy processes which have a Gaussian com- 
ponent can creep downwards. Hence for the case at hand it follows that the event 
{e Cx = x} has non-zero n-measure. 

Using the facts that W'(0+) = 2/a 2 and W(0+) = (cf. Lemmas 3.1 and 3.2) 
together with the monotonicity of W, we have that 

1 / W(x + e-e ax ) \ 

hmsup -„ (l { ,> w{x>x+e)) w{x + e) ) 

1 _ W(e) 
— nrf \ limsup»(e >x,e Gx e (x,x+e)) (3.13) 

w ( x ) eiO £ 
= 0. 

In conclusion, W"(x+) exists and 

W"(x+) = -W'(0+)n(e> Xl e ax =x)+n(e>x)W'(x), 

that is to say, 

N a 2 (W'(x) 2 „, 1 

»(£ > x, e ax = x) = — | - | . (3. 14) 

We shall now show that there exists a left second derivative W"(x—) which may 
replace the role of W"(x+) on the right hand side of (3.14) thus completing the 
proof. For this we need to recall from Doney's description of the excursion measure 
as limit (cf. [33]) that for AeJ, 

n(A,t<0=lim Py{A > t<T °\ 

yiO y 

and moreover this identity may be extended to stopping times. From this we may 
write _ 

Py(X T + =x;t+ < Tq) 
n(e > x,e a =x) = lim , (3.15) 

yio y 

where P y is the law of —X when issued from y. From part (ii) of Theorem 2.6 we 
also know that 

a 2 



P y (X z +=x;T+<~) = —W'(x-y) 



2 



42 



Alexey Kuznetsov, Andreas E. Kyprianou and Victor Rivero 



Using the Strong Markov Property, the above formula and the fact that X creeps 
upwards it is straightforward to deduce that 

P y (X z , =*;t - < r+) = x ^W'(x) 

and hence 

+ x <* 2 f i, W(x-y) ., A 

P y (X T + = x;r+ < t ") = — \^\ x -y)--^M-w'(x)y 

Returning to (3.15) we compute, 

Note that the existence of W"(x—) is guaranteed in light of (3.15). We see that the 
final equality above is identical to (3.14) but with W"(x+) replaced by W"(x—). 

Thus far we have shown that a second derivative exists everywhere. To complete 
the proof, we need to show that this second derivative is continuous. To do this, it 
suffices to show that nie > x, e 0x = x) is continuous. To this end note that a straight- 
forward computation, similar to (3.15) but making use of Part (i) of Theorem 2.7 
and L'Hopital's rule to compute the relevant limit, shows that 

n(e > x, e Cx > x) = jT j W'(x - y) - ^^W(x - y) J JT(y)dy, 

where we recall that TI(y) = fl(— °°, —y). Hence it is now easy to see that 

n(e > x, Eo x = x) = n(e > x) — n(e > x, e 0x > x) 

is continuous, thus completing the proof. □ 

Chan et al. [24] take a more analytical approach to the smoothness of scale func- 
tions by exploring its connection with the classical renewal equation. They note that, 
on account of (2.1 1), it suffices to consider smoothness in the case that y/(0+) > 
(i.e. <£(()) = 0). Within this regime, the basic idea is to start with the obvious state- 
ment that under the measure P x the Levy process will either cross downwards below 
zero by either creeping across it or jumping clear below it. Taking account of (2.29), 
(2.30) and (2.34) in its limiting form as a t °°, we thus have 
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1 



W(x) 



W<~) 

P x (X T -=0)+P x (X tQ <0) 

(J 2 y-oo 

— W(x) + {W(x)-W(x-y)}n(y)dy 



1 - 



VA'(0+) 



JO 



%-W'(x)+ [ X W'{z)n(x-z)dz 




(3.16) 



where J7(z) = n(y)dy. Hence, 



i = T w ' {x)+ Jo W ' {x - z){n{z) + VW) }dz 



(3.17) 



and after a little manipulation, we arrive at the classical renewal equation 



where f(x) = o 2 W'{x)/2, and g(x) = -2a- 2 ft TI(y)&y -2a- 2 y/' (0+)x, and mo- 
mentarily we have assumed that a 2 > 0. By engaging with the well known con- 
volution series solution to the renewal equation they establish the following result 
beyond the scope of Theorem 3.10. 

Theorem 3.11. Suppose that X has a Gaussian component and its Blumenthal- 
Getoor index belongs to [0,2), that is to say 



J\x\<\ 

Then for each q > and n = 0, 1 , 2, .. ., e C n+3 (0, °o) if and only ifU e C n (0, °o). 

In fact the method used to establish the above theorem can also be used to prove 
Theorem 3.10. It is also apparent from the analysis of [24] that their method can- 
not be applied when considering renewal equations of the form 1 = f*g, which 
would be the form of the resulting renewal equation in (3.17) when a — with an 
appropriate choice of / and g. 

However when X does have paths of bounded variation, the aforementioned 
is possible following an integration by parts and a little algebra in the convo- 
lution on the right hand side of (3.16). In that case we take / = 8W(x) and 
g(x) = 8 l /q n(y)dy where we recall that W(0+) = S~ l and 8 is the coefficient 
of the linear drift when X is written according to the decomposition (2.1). Note that 
this integration by parts is not possible if X has paths of unbounded variation with 
no Gaussian component as the aforementioned choice of g is not finite. 

The same analysis for the Gaussian case but now for the bounded variation setting 
in [24] yields the following result. 



/=!+/**, 
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Theorem 3.12. Suppose thatX has paths of bounded variation and —17 has a den- 
sity k{x), such that 7t(x) < C|x| _1_a in the neighbourhood of the origin, for some 
a < 1 and C > 0. Then for each q>0andn = 1,2,..., G C" +1 (0,°o) if and 
onlyifTl eC"(0,°o). 

Essentially the results of [24] are primarily results about the smoothness of the 
solutions to renewal (or indeed Volterra) equations which are then applied where 
possible to the particular setting of scale functions. Unfortunately this means that 
nothing has been said about the case that X has paths of unbounded variation and 
(7 = to date. Moreover, the results in the last two theorems above are subject to 
conditions which do not appear to be probabilistically natural other in so much as 
providing the technical basis with which to push through their respective analytical 
proofs. None the less they go part way to addressing Doney's conjecture for scale 
functions 1 as follows. 

Conjecture 3.13. For k = 0, 1, 2, • • • 

1. if a 2 >0then 

w e c* +3 (o,°o) e c*(o,°o), 

2. if a = and -, |x|il(dx) = °o then 

w e c* +2 (o,°o) e c*(o,°=>), 

3. if a = and /(_ 1)0 ) |jcjJT(djc) < °° then 

w e c* +1 (o,oo)«ne c k {o,°o). 

On a final note we mention that Doring and Savov [35] have obtained further 
results concerning smoothness for potential measures of subordinators which has 
implications for the smoothness of scale functions. A particular result of interest in 
their paper is understanding where smoothness breaks down when atoms are intro- 
duced into the Levy measure. 



Curious about the results in [24], R.A.Doney produced a number of specific examples of Levy 
processes whose scale functions exhibited analytical behaviour that lead to his conjecture (personal 
communication with A.E.Kyprianou). 



Chapter 4 

Engineering scale functions 



4.1 Construction through the Wiener-Hopf factorization 

Recall that the spatial Wiener-Hopf factorization can be expressed through the iden- 
tity 

y(A) = (A-#(0))*(A), (4.1) 

for A > where 0(A) is the Laplace exponent of the descending ladder height 
process. The killing rate, drift coefficient and Levy measure associated with are 
given by y/(0+) V 0, C7 2 /2 and 

/"'CO 

r(x,oc)=e*(°) JC / e-*(°)"II(-oo,-ii)dii, forx>0, 

Jx 

respectively. Moreover, from this decomposition one may derive that 

rx roo 

W(x)=e^ x / e-*^ y / dt-P(H,€dy), x>0. (4.2) 
Jo Jo 

This relationship between scale functions and potential measures of subordina- 
tors lies at the heart of the approach we shall describe in this section. Key to the 
method is the fact that one can find in the literature several subordinators for which 
the potential measure is known explicitly. Should these subordinators turn out to be 
the descending ladder height process of a spectrally negative Levy process which 
does not drift to — oo, i.e. <£(()) = 0, then this would give an exact expression for its 
scale function. Said another way, we can build scale functions using the following 
approach. 

Step 1. Choose a subordinator, say H, with Laplace exponent 0, for which one 
knows its potential measure or equivalently, in light of (2.24), one can explicitly 
invert the Laplace transform l/<j>(6). 

Step 2. Verify whether the relation 
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y(A) :=X<j>(X), X>0, 
defines the Laplace exponent of a spectrally negative Levy process. 

Of course, for this method to be useful we should first provide necessary and 
sufficient conditions for a subordinator to be the downward ladder height process of 
some spectrally negative Levy process or equivalently a verification method for the 
Step 2. 

The following theorem, taken from Hubalek and Kyprianou [48] (see also Vigon 
[104]), shows how one may identify a spectrally negative Levy process X (called the 
parent process) for a given descending ladder height process H. The proof follows 
by a straightforward manipulation of the Wiener-Hopf factorization (4.1). 

Theorem 4.1. Suppose that H is a subordinator, killed at rate K > 0, with drift 8 
and Levy measure Y which is absolutely continuous with non-increasing density. 
Suppose further that <p > is given such that q>K = 0. Then there exists a spectrally 
negative Levy process X, henceforth referred to as the 'parent process', such that 
for all x > 0, P(t+ < °°) = e~ ipj: and whose descending ladder height process is 
precisely the process H. The Levy triple (a, (7, 17) of the parent process is uniquely 
identified as follows. The Gaussian coefficient is given by <J — V28. The Levy mea- 
sure is given by 

n(-~,- x ) = pr(jc,oo) + ^(jc). (4.3) 

Finally 

a= ( xn(6x)-K (4.4) 

A— ,-i) 

if(p = and otherwise when (p > 

a =\ a2( P+yJ { _ (e^-l-^l {j[ >-i } )n(dx). (4.5) 
In all cases, the Laplace exponent of the parent process is also given by 

y/(e) = (e-<p)<Ke), (4.6) 

for9>0 where <j>(9) = -logE(e- e ^)- 

Conversely, the killing rate, drift and Levy measure of the descending ladder 
height process associated to a given spectrally negative Levy process X satisfying 
(p = <P(0) are also given by the above formulae. 

Note that when describing parent processes later on in this text, for practical 
reasons we shall prefer to specify the triple (a, 17, y/) instead of (a, <J,TI). However 
both triples provide an equivalent amount of information. It is also worth making an 
observation for later reference concerning the path variation of the process X for a 
given a descending ladder height process H. 
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Corollary 4.2. Given a killed subordinator H satisfying the conditions of the previ- 
ous Theorem, 

( i) the parent process has paths of unbounded variation if and only ifY(0, °°) = °o 
or 8 > 0, 

(ii) if Y (0,°o) = A < °° then the parent process necessarily decomposes in the 
form 

X, = (K + X-8<p)t + V28B,-S t , (4.7) 

where B = {B, :t>0}isa Brownian motion, S — {S, : t > 0} is an independent 
driftless subordinator with Levy measure V satisfying v(x,°°) = JT(— °°, —x). 

Proof. The path variation of X follows directly from (4.3) and the fact that a = 
Vl8. Also using (4.3), the Laplace exponent of the decomposition (4.7) can be 
computed as follows with the help of an integration by parts; 

(K+x-s<p)e + se 2 -<pe r e - ex r(x,oo)dx-e fV^^oOck 

Jo Jo dx 

= (K + r(0,oo) - 8<p)6 + 8G 2 -(p£(l-e- ex )^(x)dx - 6 j\- dx ^{x)Ax 

= (6-<p)^K+Se+J\l-e- ex )^(x)dxy 

This agrees with the Laplace exponent yf(d) = (9 — <p)0(0) of the parent process 
constructed in Theorem 4. 1 . □ 

Example 4.3. Consider a spectrally negative Levy process which is the parent pro- 
cess of a (killed) tempered stable process. That is to say a subordinator with Laplace 
exponent given by 

He) = K+cr(-a)((r+e) a -f), 

where a e (-1, 1) \ {0}, 7 > and c > 0. For a, J3 > 0, we will denote by 

'°* (x) = Er(^Tp)' 

the two parameter Mittag-Leffler function, see (4.22) for a precise definition of it 
and a survey of some of its properties. The following well known transforms 

J\- ex x a - l £ a . a (Xx a )dx=-^-_, (4.8) 

and 

e- 8 "A- | .v-«- | #_„ ? _„(A-'.v-«)d t = T ±— - 1, (4.9) 



Jo 



10 X-Q a 

valid for a > 0, resp. a < 0, together with the well-known rules for Laplace trans- 
forms concerning primitives and tilting allow us to quickly deduce the following 
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expressions for the scale functions associated to the parent process with Laplace 
exponent given by (4.6) such that K(p = 0. 
If0<a< lthen 

If -1 < a < 0, then 

-<px 

W(x) 



K-+cr(-a)y a 

(K + cr(-a)y a ) 2 Jo y - a ^ a \K + cr{-a)r) 

Example 4.4. Let c> 0, V > and 6 e (0, 1) and be defined by 

rv 7 r(v+A + e)' 

Elementary but tedious calculations using the Beta integral allow to prove that <j) is 
the Laplace exponent of some subordinator, H. Its characteristics are K — 0, 8 = 0, 

Y(x) := T(x,oo) = F ^ye-^+ 9 - 1 ) (e* - l) 6 " 1 , x > 0. 

Moreover, we have that fig is non-increasing and log-convex, so J7 g has a non- 
increasing density. It follows from Theorem 4. 1 that there exists an oscillating spec- 
trally negative Levy process, say X, whose Laplace exponent is y/(A) = A0(A), 
A > 0, with (7 = 0, and Levy density given by — d 2 Y /dx 2 . Using again the Beta 
integral we can obtain the potential measure of H, and as a consequence the scale 
function associated to X, which turns out to be given by 

r(v + e) o p ( r e z ('~ v ) ] 



cr(v) cr(i-e)y |/ y (e^-i)^ 6 

The integral that defines this scale function can be calculated using the hypergeo- 
metric series. The particular case where V = 1, c = T{\ + 0) appears in Chaumont, 
Kyprianou and Pardo [28] 

An interesting feature of this example is that it comes together with another ex- 
ample. Indeed, observe that the first derivative of W is given by, 



W{x) = J (eZ _ 1)1+e dz, *>0, 



which is non-increasing, convex and such that the second derivative satisfies the in- 
tegrability condition Jq(1 Ax)\W" (x)\dx < °°. So, |W"(x)| defines the Levy density 
of some subordinator. More precisely, the function defined by 
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A 



>(A):= 



Hi) 



Jo 



eT Xx dW(x) 



r(V + 0) < fd-e-^) — dx A>0 



cr(v) cr(i-e)y '(e x -i) 

is the Laplace exponent of some subordinator H* , which in turn has a Levy measure 
with a non-increasing density. Hence 

y*(X) = Xf{X) = 1 ^- y A>0, 

defines the Laplace exponent of a spectrally negative Levy process that drifts to °°. 
It can be easily verified by an integration by parts that the associated scale function 
is given by 

W*(x) = f ^- ) J o X c-^ +9 - l \^-l) 9 - 1 6z = J o X n S (z)6z, *>0. 



The method described in the previous example for generating two examples of 
scale functions simultaneously can be formalized into a general theory that applies 
to a large family of subordinators, namely that of special subordinators. 



4.2 Special and conjugate scale functions 

In this section we introduce the notion of a special Bernstein functions and special 
subordinators and use the latter to justify the existence of pairs of so called conjugate 
scale functions which have a particular analytical structure. We refer the reader to 
the lecture notes of Song and Vondracek [96], the recent book by Schilling, Song 
and Vondracek [93] and the books of Berg and Gunnar [12] and Jacobs [50] for a 
more complete account of the theory of Bernstein functions and their application in 
potential analysis. 

Recall that the class of Bernstein functions coincides precisely with the class of 
Laplace exponents of possibly killed subordinators. That is to say, a general Bern- 
stein function takes the form 

<p(9) = K + 8e+ ( (l-e- ex )r(dx), forS >0, (4.10) 

•/(o,~) 

where K > 0, 8 > and Y is a measure concentrated on (0,°°) such that f^ ^(1 A 
x)r(dx)<oc. 

Definition 4.5. Suppose that <j)(9) is a Bernstein function, then it is called a special 
Bernstein function if 
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m = ^y 0>0, (4.11) 

where <j)*(9) is another Bernstein function. In this case we will say that the Bernstein 
function <j>* is conjugate to <j). Accordingly a possibly killed subordinator is called a 
special subordinator if its Laplace exponent is a special Bernstein function. 

It is apparent from its definition that 0* is a special Bernstein function and is 
conjugate to <p*. In [46] and [97] it is shown that a sufficient condition for <p to be a 
special subordinator is that T(x,°°) is log-convex on (0,°o). 

For conjugate pairs of special Bernstein functions <j) and <j)* we shall write in 
addition to (4.11) 

$*(9) = k*+8*9 + ( (l-e- dx )r*(dx), 9>0, (4.12) 

J(0,oo) 

where necessarily Y* is a measure concentrated on (0,°°) satisfying /j 0os j(l A 
x)Y* (dx) < °o. One may express the triple (k* , 8* , Y* ) in terms of related quantities 
coming from the conjugate 0. Indeed it is known that 

(0, K>0, 

K = } (s+/ (0H xr(dx))~ ,k = o- 

5 » = f0, 5>0orr(0,oo)=oo, 
\()f + r(0,oo))- 1 ,5=0andr(0,oo)<oo. 

Which implies in particular that k*k = = 8*8. In order to describe the measure 
Y* let us denote by W(dx) the potential measure of <j). (This choice of notation will 
of course prove to be no coincidence). Then we have that W necessarily satisfies 

W(6x) = 8*8o(dx) + {K*+Y*(x,°°)}dx, forx>0, 

where 8q (dx) is the Dirac measure at zero. Naturally, if W* is the potential measure 
of <j)* then we may equally describe it in terms of (k, 8, Y). In fact it can be easily 
shown that a necessary and sufficient condition for a Bernstein function to be special 
is that its potential measure has a density on (0,°°) which is non-increasing and 
integrable in the neighborhood of the origin. 



and 



We are interested in constructing a parent process whose descending ladder 
height process is a special subordinator. The following theorem and corollary are 
now evident given the discussion in the current and previous sections. 

Theorem 4.6. For conjugate special Bernstein functions (j) and <p* satisfying (4.10) 
and (4.12) respectively, where Y is absolutely continuous with non-increasing den- 
sity, there exists a spectrally negative Levy process that does not drift to — °° whose 
Laplace exponent is described by 
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w{e) = ^W) =e ^f° r0 ^ > (4 - 14) 

and whose scale function is a concave function and is given by 

W(x) = 8* + k*x+ [ X Y*(y,oo)dy. (4.15) 
Jo 

The assumptions of the previous theorem require only that the Levy and poten- 
tial measures associated to have a non-increasing density in (0,°°), respectively; 
this condition on the potential measure is equivalent to the existence of (j>* . If in 
addition it is assumed that the potential density be a convex function, in light of the 
representation (4.15), we can interchange the roles of and 0*, respectively, in the 
previous theorem. The key issue to this additional assumption is that it ensures the 
absolute continuity of Y* with a non-increasing density and hence that we can apply 
the Theorem 4.1. We thus have the following Corollary. 

Corollary 4.7. If conjugate special Bernstein functions (j) and 0* exist satisfying 
(4.10) and (4.12) such that both Y and Y* are absolutely continuous with non- 
increasing densities, then there exist a pair of scale functions W and W*, such that 
W is concave, its first derivative is a convex function, (4.15) is satisfied, and 

W*(x) = 8 + kx+ [ X Y(y, oo)dy. (4.16) 
Jo 

Moreover, the respective parent processes are given by (4.14) and 

w * ie) = W) =orie) ' (4A7) 

It is important to mention that the converse of the latter Theorem and Corollary 
hold true but we omit a statement and proof for sake of brevity and refer the reader 
to the article [70] for further details. 

For obvious reasons we shall henceforth refer to the scale functions identified in 
(4.10) and (4.12) as special scale functions. Similarly, when W and W* exist then 
we refer to them as conjugate (special) scale functions and their respective parent 
processes are called conjugate parent processes. This conjugation can be seen by 
noting that thanks to (4.11). 

W*W*(dx) =dx. 



4.3 Tilting and parent processes drifting to — °° 

In this section we present two methods for which, given a scale function and asso- 
ciated parent process, it is possible to construct further examples of scale functions 
by appealing to two procedures. 
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The first method relies on the following known facts concerning translating the 
argument of a given Bernstein function. Let <p be a special Bernstein function with 
representation given by (4.10). Then for any j3 > the function <j>p(d) = 0(0 + 
P), 9 > 0, is also a special Bernstein function with killing term Kp = 0(J3), drift 
term dp = d and Levy measure Yp(dx) = e~P x Y(dx), x > 0, see e.g. [92] Section 
33. Its associated potential measure, Wp, has a decreasing density in (0,°°) such 
that Wp(dx) = e~P x W'(x)dx, x > 0, where W' denotes the density of the potential 
measure associated to 0. Moreover, let 0* and <j>p, denote the conjugate Bernstein 
functions of and typ , respectively. Then the following identity 

fp(e) = f (e+p)-f (l3)+l3j~(l-e- ex y-P x W'{x)dx 7 0>O, (4.18) 

holds. Note in particular that if Y has a non-increasing density then so does Yp. 
Moreover, if W' is convex (equivalently Y* has a non-increasing density) then Wp 
is convex (equivalently Yp* has a non-increasing density). These facts lead us to the 
following Lemma. 

Lemma 4.8. If conjugate special Bernstein functions and 0* exist satisfying 
(4.10) and (4.12) such that both Y and Y* are absolutely continuous with non- 
increasing densities, then there exist conjugate parent processes with Laplace expo- 
nents 

\j/p(9) = 0<j)p(0) and\j/*p(9) = 6^(6), 9 > 0, 
whose respective scale functions are given by 

Wp(x) = 8*+ J\- py T*(y,oo)dy 

= e~l }x W(x)+p [\- pz W(z)dz, 
Jo 

and 

W*p(x) = 8 + <j>(P)x + J^ y\-P z Y(dz)^j dy, (4.20) 
using obvious notation. 

All the statements in this Lemma, but the equality (4.19), follow from the previous 
discussion. The equality (4.19) is a simple consequence of the expression of W in 
(4.15) and an integration by parts. 

The second method builds on the first to construct examples of scale functions 
whose parent processes may be seen as an auxiliary parent process conditioned to 
drift to — oo. 

Suppose that is a Bernstein function such that 0(0) = 0, its associated Levy 
measure has a decreasing density and let j3 > 0. Theorem 4.1, as stated in its more 
general form in [48], says that there exists a parent process, say X, that drifts to — °° 
such that its Laplace exponent \j/ can be factorized as 



(4.19) 



x>0, 
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\i/(d) = (d-p)<i>(d), e>o. 

It follows that y/ is a convex function and y/(0) = = y/(/3), so that /3 is the 
largest positive solution to the equation y/(0) = 0. Now, let Wp be the 0-scale 
function of the spectrally negative Levy process, say Xp, with Laplace exponent 
Wp(0) := W(0 + J3)> f° r 6 > 0. It is known that the Levy process Xp is obtained 
by an exponential change of measure and can be seen as the Levy process X con- 
ditioned to drift to oo, see chapter VII in [14]. Thus the Laplace exponent \j/p can 
be factorized as \jfp(d) = 6<j>p(6) 7 for 9 > 0, where, as before, <j)p(-) := <j>(f5 + •). 
It follows from Lemma 8.4 in [67], that the 0-scale function of the process with 
Laplace exponent y/ is related to Wp by 

W(x) =?P x Wp{x), x>0. 

The above considerations thus lead to the following result which allows for the 
construction of a second parent process and associated scale function over and above 
the pair described in Theorem 4.6. 

Lemma 4.9. Suppose that (j) is a special Bernstein function satisfying (4.10) such 
that Y is absolutely continuous with non-increasing density and K — 0. Fix j3 > 0. 
Then there exists a parent process with Laplace exponent 

Y(8) = (6-PM8), 9>0, 
whose associated scale function is given by 

W(x) = 8*e Px + e Px f e^T* (y, °o)dy, x > 0, 
Jo 

where we have used our usual notation. 

In [70] the interested reader may find a discussion about the conjugated pairs 
arising in this construction. Note that the conclusion of this Lemma can also be 
derived from (2.26) and (4.15). 



4.4 Complete scale functions 

We have seen several methods that allow us to construct scale functions and pairs of 
conjugate scale functions which in principle generate large families of scale func- 
tions. In particular the method of constructing pairs of conjugate scale functions 
and its tilted versions needs the hypothesis of decreasing densities for the Levy 
measures of the underlying conjugate subordinators. This may be a serious issue 
because in order to verify that hypothesis one needs to determine explicitly both 
densities, which can be a very hard and technical task. Luckily, there is a large class 
of Bernstein functions, the so-called complete Bernstein functions, for which this 



54 



Alexey Kuznetsov, Andreas E. Kyprianou and Victor Rivero 



condition is satisfied automatically. Hence, our purpose in this section is to recall 
some of the keys facts related to this class and its consequences. 

We begin by introducing the notion of a complete Bernstein function with a view 
to constructing scale functions whose parent processes are derived from descend- 
ing ladder height processes with Laplace exponents which belong to the class of 
complete Bernstein functions. 

Definition 4.10. A function <j> is called complete Bernstein function if there exists 
an auxiliary Bernstein function 77 such that 

<p{9) = 9 2 f e-%(x)dx. (4.21) 

J(0,oo) 

It is well known that a complete Bernstein function is necessarily a special Bernstein 
function (cf. [50]) and in addition, its conjugate is also a complete Bernstein func- 
tion. Moreover, from the same reference one finds that a necessary and sufficient 
condition for <j> to be complete Bernstein is that Y satisfies for x > 

r(6x) = y^ e"^y(dy)}dx, 

where ^ ^y(dy) + ^ ^y(dy) < °°. Equivalently Y has a completely monotone 
density. Another necessary and sufficient condition is that the potential measure 
associated to <j> has a density on (0,°°) which is completely monotone, this is a 
result due to Kingman [55] and Hawkes [47]. The class of infinitely divisible laws 
and subordinators related to this type of Bernstein functions has been extensively 
studied by several authors, see e.g. [21], [103], [91], [32], [51] and the references 
therein. 

Since necessarily Y is absolutely continuous with a completely monotone den- 
sity, it follows that any subordinator whose Laplace exponent is a complete Bern- 
stein function may be used in conjunction with Corollary 4.7. The following result 
is now a straightforward application of the latter and the fact that from (4.21), any 
Bernstein function 77 has a Laplace transform §(9)/9 2 where is complete Bern- 
stein. 

Corollary 4.11. Let r\ be any Bernstein function and suppose that <j) is the complete 
Bernstein function associated with the latter via the relation (4.21 ). Write <j)* for the 
conjugate of '0 and Tj* for the Bernstein function associated with 0* via the related 
(4.21). Then 

W(x) = r)*{x)andW*(x) = r\(x), x > 0, 

are conjugate scale functions with conjugate parent processes whose Laplace expo- 
nents are given by 

^ = Jra\ = and = ^ = ^( )' 9 ^ °- 
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An interesting and useful consequence of this results is that any given Bernstein 
function 77 is a scale function whose parent process is the spectrally negative Levy 
process whose Laplace exponent is given by \j/*(d) = 9 2 /<j>(6) where is given by 
(4.21). 

Example 4.12. Let 0<a<j3<l,a,fr>0 and <f> be the Bernstein function defined 
by 

<p(e) = aeP- a +beP, e>o. 

That is, in the case where a < j3 < 1 , (j) is the Laplace exponent of a subordinator 
which is obtained as the sum of two independent stable subordinators one of param- 
eter j3 — a and the other of parameter j3 , respectively, so that the killing and drift 
term of are both equal to 0, and its Levy measure is given by 

Y(6x)=( a( fi~ a) x -(i+P-«) + ft*? x - {l+ pAdx x>0 

Hax) {r(i-p + a) x r(i-/3) ) ' >a 

In the case that j8 = 1, (j> is the Laplace exponent of a stable subordinator with 
parameter 1 — a and a linear drift. In all cases, the underlying Levy measure has a 
density which is completely monotone, and thus its potential density, or equivalently 
the density of the associated scale function W, is completely monotone. 
We recall that the two parameter Mittag-Leffler function is defined by 

^)=Lf^ W y *€R, (4.22) 

where a,j3 > 0. The latter function can be identified via a pseudo-Laplace trans- 
form. Namely, for X € R and 91(0) > A 1 /" - 7, 



/ 

Jo 



, o e e x & a p (Ax )Ax-Jq—^_ x . 



(4.23) 



The associated scale function to <j) can now be identified via 
W'(x) = ^x fS - l g a £(-ax a /b), x>0, 
which is a completely monotone function. So, the function 

\if(e) = e<p(e) = aei i - a+l +bef i+ \ e>o, 

is the Laplace exponent of a spectrally negative Levy process, the parent process. 
It oscillates and is obtained by adding two independent spectrally negative stable 
processes with stability index j3 + 1 and 1 +J3 — a, respectively. The scale function 
associated to it is given by 

W(x) = ^J X t p - 1 ^ 0l>p (-at a /b)dt, x>0. 
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The associated conjugates are given by 
and 

W*(x) = - x 1 ~P +a + - x 1_/3 x>0 (424) 

W W r(2-j3 + a) r(2-/5) ' " ( ' 

The subordinator with Laplace exponent 0* has zero killing and drift terms and 
its Levy measure is obtained by taking the derivative of the expression in (4.23). 
By Theorem 4.1 the spectrally negative Levy process with Laplace exponent \jf*, 
oscillates, has unbounded variation, has zero Gaussian terms, and its Levy measure 
is obtained by derivating twice the expression in (4.23). 

One may mention here that by letting a 1 the Continuity Theorem for Laplace 
transforms tells us that for the case 0(0) = b()P , the associated \jf is the Laplace 
exponent of a spectrally negative stable process with stability parameter 1+ j3 , and 
its scale function is given by 

The associated conjugates are given by 

0*(0)=/r 1 1 -'\ v^*(0) = fo- 1 2 - /3 , 0>O, 

and 

So that * (respectively \jf* ) corresponds to a stable subordinator of parameter 1 — j8 , 
zero killing and drift terms (respectively, to a oscillating spectrally negative stable 
Levy process with stability index 2 — J3), and so its Levy measure is given by 

n* r-oo -x) = 1H^L x p-2 x>0 

To complete this example, observe that the change of measure introduced in 
Lemma 4.8 allows us to deal with the Bernstein function 

0(0) =k{e+mf- a +b{e+mf, 0>O, 

where m > is a fixed parameter. In this case we get that there exists a spectrally 
negative Levy process whose Laplace exponent is given by 

y/(0) =kd{d + mf- a +bd{6+mf, 0>O, 



and its associated scale function is given by 
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W(x) = £ J\- m, t^- l S a ^{-at a /b)dt, x>0. 

The respective conjugates can be obtained explicitly but we omit the details given 
that the expressions found are too involved. 

We can now use the construction in Subsection 4.3. For, m,a,b>0, < a < j3 < 
1, there exists a parent process drifting to — °° and with Laplace exponent 

y/(e) = (e-m)(ae' 3 - a +/?e' 3 ), e>o. 

It follows from the previous calculations that the scale function associated to the 
parent process with Laplace exponent y/ is given by 

W & = V Jo e ~" Ut S ^i-a-t a lb)dt, x>0. 



4.5 Generating scale functions via an analytical transformation 

In Chazal et al. [69] it was shown that whenever \j/ is the Laplace exponent of a 
spectrally negative Levy process then so is 

where y/W(w) = y(u) - q, 5,j3 > and \jfW(0+) = q = if J3 = 0. Note that 
is the Laplace exponent of a spectrally negative Levy process possibly killed at 
an independent and exponentially distributed random time with rate q. In the usual 
way, when q = we understand there to be no killing. Moreover, the characteris- 
tics of the Levy process with Laplace exponent y/ are also described in the 
aforementioned paper through the following result. 

Proposition 4.13. Fix 8, j3 > with the additional constraint that y/ <? )'(0+) = q = 
iffi = 0. If has Gaussian coefficient <J and jump measure FI then ^sfi a ^ so 
has Gaussian coefficient a and its Levy measure is given by 

^ x n{Ax) + 8^ x Tl{x)Ax + 8^ x dx on (—,0), 
where TI{x) = TI{— °°, —x). 

In particular we note that ^s.pY^ is the Laplace exponent of a spectrally negative 
Levy process without killing, i.e. ^s.pW^H^) = 0. 

It turns out that this transformation can be used in a very straightforward way 
in combination with the definition of g-scale functions to generate new examples. 
Note for example that for a given y/ and ^>0we have 
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jTe-*W(«>(*)dx = for 9 > <P(q), 

and hence it is natural to use this transformation to help find the Laplace inverse (if 
it exists ) of 

1 



for j3 sufficiently large to give an expression for sV («)> the 0-scale function as- 
sociated with the spectrally negative Levy process whose exponent is ^p t sW^- The 
following result, taken from [69] does exactly this. Note that the first conclusion in 
the theorem gives a similar result to the conclusion of Lemma 4.8 without the need 
for the descending ladder height to be special. 

Theorem 4.14. Let x, J3 > such that y/(0+) = q = iffi = 0. Then, 

W ^ P M X ) = e-^W^W+jS J\-Py\V^(y)dy. (4.25) 
Moreover, if\j/(0+) < 0, then for any x, 8,q > we have 

Proof. The first assertion is proved by observing that 

\- ex w a 



L 



.^./sv w,v > 0y/(0 + j3) 



1 + P 
V/(e + j3) + 0W0 + i3) : 



which agrees with the Laplace transform of the right hand side of (4.25) for which 
an integration by parts is necessary. As scale functions are right continuous, the 
result follows by the uniqueness of Laplace transforms. 

For the second claim, first note that ^ >e y/^ = (e + <£(O)-<5)y/(0 + <i>(O))/(e + 
0(0)). A straightforward calculation shows that for all 9 + 8 > 0(0), we have 

r e -fl* e (*(o)-«w (x)dx = d + 8 

The result now follows from the first part of the theorem. □ 

Below we give an example of how this theory can be easily applied to generate 
new scale functions from those of (tempered) scale functions. 

Example 4.15. Let 

^(9) = {9+c) a -c a -qfov9> -c, 
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where 1 < a < 2 and q,c > 0. This is the Laplace exponent of an unbounded vari- 
ation tempered stable spectrally negative Levy process £, killed at an independent 
and exponentially distributed time with rate q. In the case that c — 0, the underlying 
Levy process is a spectrally negative a-stable Levy process. In that case it is known 
that 

f°° 1 

/ z- e *x a - l g a . a {qx a )Ax=- , 

Jo a — q 

and hence the scale function is given by 

W {q) {x)=x a - l g a , a (qx a ), 

for x > 0. (Note in particular that when q = the expression for the scale function 
simplifies to r(a)~ l x a ~ l ). Since 

poo 1 

L ^ cx \r^ x)Ax = ie + cr-c«-q ' 

it follows that 

W (q) (x) = e"°W r 9+c «, (x) = e- cx x a -^ a<a ((q + c a )x a ). 
Appealing to the first part of Theorem 4.14 we now know that for j3 > 0, 

W < q) (x)=e-^ x x<*- 1 £ a , a ((q+c a )x a )+l3 [\-^y y a - l g a . a ((q+c a )y a )&y. 

J $ Jo 

Note that ^^'(0+) — ac a ~ l which is zero if and only if c = 0. We may use the 
second and third part of Theorem 4.14 in this case. Hence, for any 8 > 0, the scale 
function of the spectrally negative Levy process with Laplace exponent ^gfiV^ is 

where we have used the recurrence relation for the Gamma function and r(a,b) 
stands for the incomplete Gamma function of parameters a,b > 0. Moreover, we 
have, for any j3 > 0, 

^vr w - nh) (^ r («-'-'")-^"^ r (a-i,&)) ■ 

Finally, the scale function of the spectrally negative Levy process with Laplace ex- 
ponent ^8 Yq^ is given by 
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where we have used the notation 

r (x; an + p)q" 



„to na» + £) 



4.6 Shifted scale functions 

In this final section, we make some remarks concerning how, in certain circum- 
stances, simply adding unity to an existing scale function creates a new scale func- 
tion of a completely unrelated spectrally negative Levy process. 

Theorem 4.16. Suppose that W is a complete scale function belonging to a spec- 
trally negative Levy process of unbounded variation, and whose Laplace exponent 
y/ satisfies y/(0+) = 0. Then 1+W(x) is a scale function of a spectrally negative 
Levy process, say X s , whose Laplace exponent is given by 

V{ e):=-^-for9>0. 

Moreover, it is necessarily the case that X s has paths of bounded variation and its 
jump measure, II s , satisfies 

n s (-o° 7 - x ) = w'(x) 

forx > where W is the scale function associated with the spectrally negative Levy 
process whose Laplace exponent is given by yf(6) + 6. 

Proof. The assumptions in the theorem allow us to write \jf(d) = d(j)(d) where is 
a complete Bernstein function satisfying <j> (0) = 0. 

Let W be the scale function of the spectrally negative Levy process with Laplace 
exponent + \j/(d). The key to the proof is to consider a compound Poisson sub- 
ordinator, say S — {S t : t > 0} with unit arrival rate and jump distribution whose 
jump measure is given by W(6x). Note that the latter is a probability distribution on 
account of the fact that for x > 

W[0,oo)=lim/ e~ Px W(6x)=lim—Jr — - = 1. 
1 ' Wo ./[oh wo 0(0) + 1 



Suppose that V (dx) is the renewal measure associated with this compound Poisson 
subordinator. In particular, a straightforward computation shows that 
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V(dbc) = f dt-P(S t edx) 
Jo 

Jo to » ! 
= fw*"(dx), 

n=0 

where we understand W*°(dx) = <5o(dx). 

Taking Laplace transforms we find, for j3 > 0, that 

However, since 

/ e -^{5o(dx)+W(dx)} = l + -^, 



.+-1 



we deduce that 

V(x) = l+W(x). 

Our objective is now to show that V is the scale function of the process X s . To 
this end, note that the compound Poisson subordinator behind the measure V has 
Laplace exponent 



poo 



1 - e- e *)W(dx) = 1 - -^tttt = T^TTT- 
' K ' 1+0(0) 1 + 0(0) 



Since the jump measure W is also the renewal measure associated with the Bernstein 
function 1 + 0, which is a complete Bernstein function (because is), it follows 
that W' is non-increasing and we may apply Theorem 4.1 to conclude that there 
exists a spectrally negative Levy process whose descending ladder height is the 
aforementioned compound Poisson subordinator. Moreover, this spectrally negative 
Levy process, X s , has Laplace exponent 

and its scale function is given by V = 1 + W. 

Note that the process X s necessarily has paths of bounded variation as its scale 
function has a discontinuity at the origin. □ 

Example 4.17. Consider the case of a spectrally negative stable process with index 
a € (1,2). Its Laplace exponent is given by y/(0) = 9 a and its scale function is 
given by W(x) = x a ~ l /r(a). The above theorem predicts that 1 +x a ~ 1 /r(a) is a 
scale function of a process whose Laplace exponent is given by 

al+a q2 

r{ 8)= 6 6 



e + e a e 2 - a + e' 
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Moreover, since it is known from Furrer [42] that the scale function associated with 
spectrally negative Levy process whose Laplace exponent is given 9 + 9 a is given 
by 

W(x) = l-g a _u(-x a - 1 ), 

this gives us the associated jump measure of X s . Note that this example can also be 
recovered from (4.24) with an appropriate choice of constants. 



Chapter 5 

Numerical analysis of scale functions 



5.1 Introduction 

The methods presented in the previous chapters for producing closed form expres- 
sions for scale functions are generous in the number of examples they provide, but 
also have their limitations. This is not surprising, since the scale function is defined 
via its Laplace transform, and in most cases it is not possible to find an explicit ex- 
pression for the inverse of a Laplace transform. Our main objective in this chapter is 
to present several numerical methods which allows one to compute scale functions 
for a general spectrally negative Levy process. As we will see, these computations 
can be done quite easily and efficiently, however there are a few tricks that one 
should be aware of. Our second goal is to discuss two very special families of Levy 
processes, i.e. processes with jumps of rational transform and Meromorphic pro- 
cesses, for which the scale function can be computed essentially in closed form. 

The problem of numerical evaluation of the scale function and other related quan- 
tities has received some attention in the literature. In particular, Rogers [90] com- 
putes the distribution of the first passage time for spectrally negative Levy processes 
by inverting the two-dimensional Laplace transform. The main tool is the discretiza- 
tion of the Bromwich integral and the application of Euler summation in order to im- 
prove convergence. We will describe the one dimensional version of this method in 
Section 5.3.2. Surya [100] presents an algoritm for evaluating the scale function us- 
ing exponential dampening followed by Laplace inversion; the latter performed in a 
similar way as Rogers [90]. In a recent paper Veillette and Taqqu [105] compute the 
distribution of the first passage time for subordinators using two techniques. These 
are the discretization of the Bromwich integral and Post-Widder formula coupled 
with Richardson extrapolation. Albrecher, Florin and Kortschak [5] develop algo- 
rithms to compute ruin probabilities for completely monotone claim distributions, 
which is equivalent to computing the scale function W(x) due to relation (2.8). 

In this section we will present some general ideas related to numerical evaluation 
of the scale function. We adopt the same notation as in the previous chapters. How- 
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ever, in order to avoid excessive technical details, we shall impose the following 
condition throughout. 

Assumption 1. The Levy measure fl has at most a finite number of atoms when X 
has paths of bounded variation. 

Recall that the scale function (x) is defined by the Laplace transform identity 



We know from Lemma 2.4 and Corollary 2.5 that W^'(x) exists and is continuous 
everywhere except when X has paths of bounded variation, in which case the deriva- 
tive does not exist at any point x such that an atom of of the Levy measure occurs at 
—x. Since we have assumed that there exists just a finite number of these points, we 
can apply standard results (such as Theorem 2.2 in [30]) and conclude that (x) 
can be expressed via the Bromwich integral 



where c is an arbitrary constant satisfying c> 0(q). 

In principle one could use (5.2) as a starting point for numerical Laplace inver- 
sion to give W^(x). However, this would not be a good approach from the nu- 
merical point of view. The problem here is the exponential factor e cx , which can 
be very large and would amplify the errors present in the numerical evaluation of 
the integral. Ideally we would like to choose c to be a small positive number, but 
this is not possible due to the restriction c > <l>{q). Therefore this method would 
be reasonable from the numerical point of view only when q = and <£(()) = 0, 
in which case the function W^ q \x) does not increase exponentially fast. Indeed in 
such cases there is at most linear growth. This follows on account of the fact that 
when q — <P(0) = the underlying Levy process does not drift to — °° and W is the 
renewal measure of the descending ladder height process; recall the discussion pre- 
ceding (2.25). Thanks to (2.25), this in turn implies that W(x) is a renewal function 
and hence grows at most linearly as x tends to infinity. 

In all other cases we would have to modify (5.2) in order to remove the exponen- 
tial growth of W ( ? ) (x) . It turns out that this can be done quite easily with the help of 
the density of the potential measure of the dual process X = —X, defined as 




(5.1) 




(5.2) 




(5.3) 



We know from Theorem 2.7 (iv), this function satisfies 



W {q) {x) 



v' (<?(<?)) 



-u 



u iq) (x), x>0 



(5.4) 



The theory of scale functions for spectrally negative Levy processes 



65 



therefore u^ q \x) is continuous on (0,°°) and has finite left and right derivatives at 
every point x > 0. Theorems 3.10, 3.11 and 3.12 give us more information on the 
relation between the Levy measure of X and the smoothness properties of u^ q \x). 
As we see from identity (5.4), the problem of computing the scale function (x) 
is equivalent to that of computing u^ q \x). Dealing with the latter turns out to be an 
easier problem from the numerical point of view, providing that q > and <&{q) > 0, 
since in this case is bounded. To see why this is the case, note from the earlier 
representation of in (3.5) together with (2.3), we have that 

e <P(q)x 

v'(<£(<?)) 

*(,), — i P? {q) {X„>0) 



< 



4— Pf (9) (X„<0) 

p <P(q)x 

— pf (9) (Tn <-) 

V'^fe)) v {T ° <oo} 
1 

V(*(?))' 



thereby showing boundedness. Note that similar computations to the above can be 
found in Takacs [101] and Bingham [20]. 

We need to characterize the Laplace transform of u^ q \x). The proof of the next 
proposition follows easily from (5.4). 

Proposition 5.1. Assume that <P(q) > 0. Then for Re(z) > 

CO 

Je- zx u^{x)dx = F^{z), (5.5) 
o 

where 

Fiq) {Z) = Xlf(<P(q))(z-<P(q)) WTq ' ^ 
Let us summarize our approach to computing the scale functions. When q = and 
0(0) = we will work with the scale function itself and will use (5.2) as the starting 
point for our computations. As earlier noted, in this case W { - q \x)=W{x) grows at 
worst linearly fast as x — > +°° and we do not need to worry about amplifications of 
numerical errors. When q > and <&(q) > the scale function grows exponentially 
fast as x — > +°°, thus it is better to work with the density of the potential measure 
u" q \x), and then to recover the scale function via relation (5.4). For convenience 
however we shall restrict ourselves to the latter case with the following blanket 
assumption. 
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Assumption 2. For q>0,<P(q)> 0. 



The analysis for the case q = <P(q) = 0, where we work directly with the scale 
function, is no different. 

As we will see later, our algorithms will require the evaluation of F^ q \z) for 
values of z in the half-plane Re(z) > 0. From (5.6) we find that F^(z) has a remov- 
able singularity at z = <£(<?)• It is not advisable to compute F^ q \z) via (5.6) when 
z is close to <£(<7), as this procedure would involve subtracting two large numbers, 
which would cause the loss of accuracy. There are two solutions to this problem. 
One should either make sure that z is never too close to <P(q) or alternatively one 
should use the following asymptotic expression for F^ q '(z). 

Proposition 5.2. Define a„ = y/ n ) (<t>{q)) where y/ n ) is the n-th derivative of y/. 
Then as z ^ &(q) 



I a\ 



1 A3 \a\ 



(z-<P(q)) + 0((z-<P(q)) 2 ). (5.7) 



Proof. The proof follows easily by writing down the Taylor expansion of the right- 
hand side in (5.6) centered at z = <P(q). 

As we have seen in Section 2.5 and in many other instances, virtually all fluctu- 
ation identities for spectrally negative Levy processes can be expressed in terms of 
the following three objects: the scale function (x), its derivative W^'(x) and a 
function (x) defined by (2.27), which is essentially an indefinite integral of the 
scale function. Therefore, in addition to the scale function itself, it is also important 
to be able to compute its derivative and indefinite integral. It turns out that computa- 
tion of all these quantities can be done in exactly the same way. Again, the first step 
is to remove the exponential growth as x — » +°° and express everything in terms of 
the potential density (x) as follows 

W{q) '^=^w- a(9) ' {x) ' (5 ' 8) 

and 

where we have defined 

v («)(jc)= f'aM(y)dy. (5.10) 
Jo 

We see that the problem of computing W^'(x) and Z( q \x) is equivalent to the 
problem of computing u^'(x) and v^)(x). The following result is an analogue of 
proposition 5.1 and it gives us Laplace transforms of «W'(x) and v^ q \x). 

Proposition 5.3. Assume that q>0 and <P(q) > 0. Then for Re(z) > 



u 



(<?)' 



(x)e--dx = zF<«) (z) + (0+) - — l — - , (5.11) 



o W w V ' ¥(<P(q)) 
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and 



(5.12) 



z 



Proof. The proof follows from proposition 5.1 and integration by parts. 

Now we have reduced all three problems to an equivalent "standardized" form. 
Equations (5.4), (5.8), (5.9) and Propositions (5.1) and (5.3) show that the problem 
of computing of W^(x), W^'(x) and Z^ q \x) is equivalent to a standard Laplace 
inversion problem. One has an explicit expression for the function g(z), which is 
analytic in the half -plane Re(z) > and is equal to the Laplace transform of f(x) 



and one wants to compute the function f(x) for several values of x > 0. From now 
on we will concentrate on solving this problem. 

This Laplace inversion problem has been well-studied and it has generated an 
enormous amount of literature. We will just mention here an excellent textbook by 
Cohen [30], very helpful reviews and research articles by Abate, Choudhurry, Whitt 
and Valko [4, 1, 2] and works of Filon [39], Bailey and Swarztrauber [11] and Iserles 



This chapter is organized as follows. In Sections 5.2 and 5.3 we present four 
general numerical methods for computing f(x). The first approach starts with the 
Bromwich integral and is based on Filon's method coupled with fractional discrete 
Fourier transform and Fast Fourier Transform techniques. The next three methods, 
i.e. Gaver-Stehfest, Euler and Talbot algorithms, also start with the Bromwich inte- 
gral, but the discretization of this integral is done in a different way, and in every 
case (except for Talbot method) some acceleration procedure is applied. These three 
methods typically require multi-precision arithmetic. In Sections 5.4 and 5.5 we 
present two families of processes for which the scale function can be computed ex- 
plicitly, as a finite sum or an infinite series, which involve exponential functions, 
derivative of the Laplace exponent y/(z) and the solutions to equation y/(z) = q. In 
this case computing the scale function can be done in an extremely efficient and 
accurate way, and later we will use these processes as benchmarks to test the perfor- 
mance of the four general methods. In Section 5.6 we discuss the results of several 
numerical experiments and in Section 5.7 we present our conclusions and provide 
some recommendations on how to choose the right numerical algorithm. 



5.2 Filon's method and fractional fast Fourier transform 




(5.13) 



[49]. 



The starting point for Filon's method is an expression for f(x) which identifies it as 
a Bromwich integral. That is, 
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/W = ^i / sto^dz, *>0 



(5.14) 



c+iR 



where c is an arbitrary positive constant. Since we are only interested in f(x) for 
positive values of x, e equation (5.14) can be written in terms of the cosine transform 
as follows, 



Now our plan is to compute the above integral numerically. 

The integral in (5.15) is an oscillatory integral, which means that the integrand 
is an oscillating function. Evaluating oscillatory integrals is not as straightforward 
as it may seem, and one should be careful to choose the right numerical method. 
For example, in (5.15), the cosine function can cause problems. When x is large 
it oscillates rapidly with the period 2n/x. Therefore if we simply discretize this 
integral using the trapezoid or Simpson's rule, we have to make sure that the spacing 
of the discretization h satisfies h<sc.2%jx. When x is large, this restriction forces us to 
take h to be a very small number, therefore we need a huge number of discretization 
points and the algorithm becomes slow and inefficient. The main benefit of Filon's 
method is that it helps one to avoid all these problems. 

Filon's method applies to computing oscillatory integrals over a finite interval 
[a,b] of the form 



Here & C G stands for cosine transform of the function G. In order to describe the 
intuition behind Filon's method, let us revisit Simpson's rule. It is well known that 
Simpson's rule for the integral in (5.16) can be obtained in the following way. We 
discretize the interval [a,b], approximate function u G(u) cos(ux) by the second 
order Lagrange interpolating polynomials in the M-variable on each subinterval, and 
then integrate this approximation over [a,b\. Filon's method goes along the same 
steps, except for one crucial difference. The function G(u) is approximated by sec- 
ond order Lagrange interpolating polynomials, which is then multiplied by cos(wx) 
and integrated over [a,b\. Due to the fact that the product of trigonometric functions 
and polynomials can be integrated explicitly, we still have an explicit formula. In do- 
ing this we separate the effects of oscillation and approximation. The approximation 
for the function G(u) is usually quite smooth and does not change very fast. Then 
any effect of oscillation disappears since we compute the integral against cos(mx) 
explicitly. 

Let us describe this algorithm in full detail. Our main references are [39], [41] 
and [49]. We take N to be an integer number and define h = (b — a)/(2N) and 
u„ = a + nh, < n < 2N. Then we denote g n — G(u n ), < n < 2N and introduce 
the vectors u = [uq,u\,.. . ,U2n] an d g = [go,gi: ■ ■ ■ igiN} - For fee {1,2} we define 




(5.15) 




(5.16) 



AT-1 

%(g,»,x) = £ g2n+kCOs(xu 2 „ +k )- 
n=0 



(5.17) 
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Next, on each subinterval [«2«,«2n+2], < n < N we approximate G(u) by a La- 
grange polynomial of degree two, which gives us a composite approximation of the 
following form 



N-i 



G(u;N) = £ 1{h 2 „<„<« 2 „ +2 } 

n=0 

1 



g2n+\ + (g2n+2 ~ g2n) (« ~ «2«+l ) 



+ ^2-fe«+2-2§2n+l +g2n)(M-«2«+l) 2 



Multiplying this function by cos(mx), integrating over the interval [a,b] and simpli- 
fying the resulting expression we obtain the final expression for Filon's method, 



V 

& c G{x;N) = J G(u;N)cos(ux)du 



(5.18) 



hA(hx)(G(b) sin(fax) — G(a) sin(ax)) 



+ hB(hx) 
+ hC{hx)^i{g,u,x), 



^2(g,u,*) — - (G(b)cos(bx) — G(fl)cos(ax)) 



where 



A(0) 
B(9) =2 
C(0) =4 



1 sin(20) 2sin(0) 2 



e 2e 2 e 3 ' 

l+cos(0) 2 sin(20) 



2 3 

sin(0) cos(0)~ 



3 



02 



(5.19) 
(5.20) 
(5.21) 



Note that by construction, Filon's approximation is exact for polynomials of de- 
gree two or less. It is also known that the error of Filon's approximation is 0(h 3 ), 
provided that G^(u) is continuous, see [41]. 

We see that in order to evaluate Filon's approximation (5.18) we have to compute 
two finite sums ^(g,u,x) defined by (5.17). If we want to compute & c G{x;N) 
for just a single value of x then this obviously requires 0(N) operations. By the 
same reasoning, if we want to compute JF c G(x;N) for N equally spaced points x = 
[xo,x\,...,xn-i], then we would have to perform 0(N 2 ) operations. However, the 
special structure of these sums makes it possible to perform the latter computations 
in just 0(Nln(N)) computations. The main idea is to use fast Fourier transform 
(FFT). This works as follows. Assume that we want to compute J^ c G(x;N) forN x = 
N values x m = xo + mS x , < m < N. We rewrite the finite sums in (5.18) in the 
following form 
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%(g,U,x m ) =Re 



N-l 

\(a+kh)x m ^ g k ne ,i(2fc5j)nm 



K=0 



(5.22) 



where we have defined g^ „ = g2n+k e>2hnx ° ■ Then, assuming that parameters h and 5* 
satisfy 

h8 x = ^, (5.23) 

the expression in the right-hand side of (5.22) is exactly in the form of the discrete 
Fourier transform (see [11]), and it is well-known that it can be evaluated for all 
m = 0, l,..,N— 1 in just 0(N\n(N j) using the fast Fourier transform technique. 

The restriction (5.23) is quite unpleasant as it does not allow us to choose the 
spacing between the discretization points in the M-domain in (5.16) independently 
of the spacing in the x-domain. However there is an easy solution to this problem, 
namely the fractional discrete Fourier transform (see [11]), which is defined as a 
linear transformation, which maps a vector v = [vo, Vi, . . . , v^-i] into a vector V = 

\y ,v u ...,v N -i] 

N-l 

V m =I> n e i0! ™ m = 0,l,...,AT-l. (5.24) 

n=0 

It turns out that for any ael one can still compute the values of V m for m = 
0,1,.., AT- 1 in just 0(N\n(N)) operations, see [11] for all the details. 

Let us summarize the main steps of algorithm. First, choose a small number 
c > 0, so that the factor e cx is not too large for the values of x that interest us. Then 
set a — and choose the value of the cutoff, i.e. a large number b > such that the 
integral 

/■CO 

/ Re [g(c + it/)]cos(t«)di< 
Jb 

is sufficiently small. Then, choose the number of discretization points N in the u- 
domain and define u n — nb/(2N), < n < 2N and g n = Re[g(c + iu n )]. Choose 
8 X and xq and compute the approximation (5.18) using the fractional Fast Fourier 
Transform. This gives us N values of /(xq + m8 x ) for < m < N, with the error 
bound 0(N~ 3 ), at the computational cost of 0(N\n(N)) operations. 

There is another trick that might be very useful when implementing Filon's 
method. In many examples the integrand Re[g(c + iM)] in (5.15) changes quite 
rapidly when u is small while it changes slowly for large values of u. This means 
that we would have better precision and would need fewer discretization points if 
we were able to place more of them near u = and fewer of them for large u. This 
can be easily achieved by dividing the domain of integration 

I Re{g(c + iu)]cos(ux)du= Y\ / Re[g(c + iM)]cos(t«)dM, (5.25) 
Jo PoJbj 
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where = bo < b\ < ■■■ < b n = b. Each integral over [bj,bj + i] can be evaluated 
using Filon's method to produce results on the same grid of the x-variable. Choosing 
bj so that the spacing bj + \ —bj increases allows us to concentrate more points where 
they are needed (near u = 0) and fewer points in the regions far away from u = 0. 



5.3 Methods requiring multi-precision arithmetic 

The methods presented in this section can give excellent performance, but the price 
that one has to pay is that they all require multi-precision arithmetic. Our main 
references for this section are [4], [2], [1] and [30]. These methods are grouped 
together since they all give a similar expression for the approximation to f{x) 



where the coefficients a n and b n depend only on M and do not depend on functions 
/ and g. 



5.3.1 The Gaver-Stehfest algorithm 

The first method that we discuss is the Gaver-Stehfest algorithm, see [4] and Section 
7.2 in [30]. This algorithm is based on the Gaver's approximation [43], which can 
be considered as a discrete analogue of the Post-Widder formula (see Section 2.3 in 



It turns out that both Gaver's and Post-Widder formulas have a very slow conver- 
gence rate, therefore one has to apply some acceleration algorithm, such as Salzer 
transformation for the Gaver's formula, which was proposed by Stehfest [98], or 
Richardson extrapolation for Post-Widder formula which was used by Veillette and 
Taqqu [105]. 

We refer to [4] for all the details and background on the Gaver-Stehfest algorithm, 
here we just present the final expression. The function f(x) is approximated by 
f GS (x;M), which depends on a single integer parameter M and is defined as 




[30]) 




ln(2) 



2M 



f GS (x;M) = 



x 



Y^a n g(n\n(2)x l ) , 



(5.26) 



n=\ 



where the coefficients a n are given by 
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Note that the coefficients a n are defined as finite sums of possibly very large num- 
bers, and that these coefficients have alternating signs. This means that we will 
lose accuracy in (5.26) due to subtracting very large numbers, thus we have to use 
multi-precision arithmetic to avoid this problem. Abate and Valko [2] (see also [4]) 
recommend the following "rule of thumb". If we want j significant digits in our 
approximation, we should set M = [1.1/1 (the least integer greater than or equal to 
1.1 7) and set the system precision at [2.2M] . It is also useful to check the accuracy 
of computation of a„ using the fact that 

2M 

j>„=o. 

n=l 

We see that Gaver-Stehfest algorithm should have efficiency around 0.9/2.2 <~ 
0.4, which is the ratio of the system precision to the number of significant digits 
produced by the approximation. While algorithms presented below all have higher 
efficiency, the Gaver-Stehfest algorithm has a big advantage that it does not require 
the use of complex numbers. This can improve performance, since it is faster to 
perform computations with real numbers than with complex numbers. This feature 
is also useful, because not every multi-precision software has a built-in support for 
complex numbers. 



5.3.2 The Euler algorithm 



The next two algorithms are based on the Bromwich integral representation, which 
follows from (5.14) after the change of variables of integration u ^ v/x 



/(*) = 



2nx 



e lv dv, xe 



(5.27) 



Note that we have also changed c i->- c/x. The main idea behind Euler algorithm 
is to approximate the integral in (5.27) by a trapezoid rule and then apply Euler 
acceleration method to improve the convergence rate. Again, all the details can be 
found in [1] and [4]. We present here only the final form of this approximation 



f E (x;M)= 1 ^£(-l)"fl„Re[g((ln(l0")+ra«)x- 1 ) 



(5.28) 



where the coefficients a n are defined as 
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Note that while coefficients a n are real, formula (5.28) still requires evaluation of 
g(z) for complex values of z. 

In the case of Euler algorithm, Abate and Whitt [4] recommend to set M = [1.7 j] 
if j significant digits of are required, and then set the system precision at M. Coeffi- 
cients a„ can be precomputed, and the accuracy can be verified via condition 



We see that the efficiency of the Euler algorithm should be around 1/1.7 ~ 0.6. 



5.3.3 The fixed Talbot algorithm 

The Talbot algorithm [102] also starts with integral representation (5.27), but then 
the contour of integration is transformed so that Re(z) — > — °° on this contour. Note 
that Re(z) is constant on the contour of integration in the Bromwich integral (5.27). 
This transformation of the contour of integration has a great benefit in that the inte- 
grand g(z) exp(zx) converges to zero much faster. On the negative side, this method 
only works when (i) g(z) can be analytically continued into the domain |arg(z) | < n, 
and (ii) g(z) has no singularities far away from the negative half-line. In our case 
these two conditions are satisfied for processes whose jumps have completely mono- 
tone density, as in this case all of the singularities of (z) lie on the negative half- 
line. Meromorphic Levy processes, presented in Section 5.5, exhibit this property 
for example. It is also worthy of note for future reference that such processes are 
dense in the class of processes with completely monotone jumps. As we will see 
later, Talbot method gives an excellent performance for processes with completely 
monotone jumps. 

Again, we refer to [4], [2] and [30] for all the details, and present here only the 
final form of the approximation 



2M 



L(-i)»fl»=o. 



«=0 



f{x;M) = - £ lm[a n g(b n x- 1 )] 



j M-l 



(5.29) 



X n=0 



where 




and 
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For this algorithm, Abate and Valko [2] (see also [4]) recommend using the same 
precision parameters as for the Euler algorithm: one should set M = [1.7 j~\ if j sig- 
nificant digits of are required, and then set the system precision at M. The efficiency 
of this algorithm should also be close to 0.6. 



5.4 Processes with jumps of rational transform 

It is well known that the Wiener-Hopf factorization and many related fluctuation 
identities can be obtained in closed form for processes with jumps of rational trans- 
form. The Wiener-Hopf factorization for the two-sided processes having positive 
and/or negative phase-type jumps (a subclass of jumps of rational transform) was 
studied by Mordecki [81], Asmussen, Avram and Pistorius [7] and Pistorius [85], 
while Wiener-Hopf factorization for a more general class of processes having posi- 
tive jumps of rational transform was obtained by Lewis and Mordecki [80]. Expres- 
sions for the scale functions for spectrally-negative processes with jumps of rational 
transform follow implicitly from these papers, their explicit form was obtained in a 
recent paper by Egami and Yamazaki [37]. 

In order to specify these processes, let us define the density of the Levy measure 
as follows 

m 

n(x) = l {x<0] ^aj\xp- l e PJ x , (5.30) 
y=i 

where rrij e N and Re(p 7 ) > 0. It is easy to see that the Laplace exponent is given 
by 

a 1 m r i 

Y(z) = -z 2 + mz+I>K-1)! (pj-z) "y- P j" ,J (5.31) 

and that is a rational function 

(5.32) 

where deg(g) = M = £? =1 mj and deg(P) = N, where N = M + 2 (resp. M + 1) 
if a > (resp. a — 0). The next proposition gives us valuable information about 
solutions of the equation \j/(z) = q. 

Proposition 5.4. 

(i) For q > or q = and y/(0) < equation \)f(z) = q has one solution 
z = in the half-plane Re(z) > and N — 1 solutions in the half-plane 
Re(z) < 0. 

(iij For q = and y/(0) > (resp. y/(0) =0) equation = q has a solution 
Z = of multiplicity one (resp. two) and N — 1 (resp. N — 2) solutions in the 
half-plane Re(z) < 0. 
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( iii) There exist at most M + N — 1 complex numbers q such that the equation 
\jf(z) — q has solutions of multiplicity greater than one. 

Proof. The proof of (i) and (ii) is trivial and follows from (5.32) and the general 
theory of scale functions. Let us prove (iii). Assume that the equation y/(z) = q 
has a solution z = Zo of multiplicity greater than one, therefore y/(zo) = 0. Using 
(5.32) we find that y/(z ) = implies P'{z )Q(zo) - P(z Q )Q(zo) = 0. The polyno- 
mial H(z) = P'{z)Q{z) — P{z)Q{z) has degree M + N—l, thus there exist at most 
M + N — 1 distinct points Zk for which yf'(zk) = 0, which implies that there exist at 
most M+N — 1 values q, given by qk = \j/(zk), for which the equation \jf(z) = q has 
solutions of multiplicity greater than one. 

The statement in Proposition 5.4 (iii) is quite important for numerical calcula- 
tions. While it's proof is quite elementary, we were not able to locate this result 
in the existing literature. The implication of this result is that for a generic Levy 
measure defined by (5.30) and general q > it is extremely unlikely that equation 
\jf(z) = q has solutions of multiplicity greater than one. So for all practical purposes 
(unless we are dealing with the case y/(0) = q = 0) we can assume that all the 
solutions of y/(z) = q have multiplicity equal to one, as we will see later this will 
considerably simplify all the formulas. 

Computing the scale function, or equivalently, the potential density u^ (x) is a 
trivial task for this class of processes. Let us consider the general case, and assume 
that equation \jf(z) = q has n distinct solutions — £i, — £2, — in the half -plane 
Re(z) < 0, and the multiplicity of z = is equal to nj. Due to (5.32) it is clear that 
N —1 = YIj=i n j- Rewriting the rational function 1 / (y/(z) — q) as partial fractions 

1 1 » 1 » 1 C j 

W^q = V '(0(q))(z-0(q)) + £ £ JzTW ' ^ 

and using Proposition 5. 1 we can identify u^ (x) as follows 

^(x^-fe-^jr^^- 1 , *>0. (5.34) 

i=\ k=\ y K L >- 

Note, that if all the solutions have multiplicity one, i.e. nj = 1 and n + 1 = N = 
deg(P), then (5.33) implies that 

1 

and we have a much simpler expression for the potential density 

N-l -Qx 

u^(x) = -^-—- x>0. (5.35) 

,=1 Y \ hj) 
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We would like to stress again that Proposition 5.4 (iii) tells us that fi( q \x) can be 
computed with the help of (5.35) for all but a finite number of q, unless we are 
dealing with the case y/(0) = q = 0. 



5.5 Meromorphic Levy processes 

The main advantage of processes with jumps of rational transform is that the nu- 
merical computations are very simple and straightforward. Everything boils down 
to solving a polynomial equation \j/(z) = q and performing a partial fraction decom- 
position of a rational function l/(y/(z) - q). On the negative side, it is clear that 
these processes can only have compound Poisson jumps, while in applications it is 
often necessary to have processes with jumps of infinite activity or even of infinite 
variation. Meromorphic Levy processes, which were recently introduced in [65], 
solve precisely this problem. They allow for much more flexible modeling of the 
small jump behavior, yet all the computations can be done with the same efficiency 
as for processes with jumps of rational transform. 

In order to define a spectrally-negative Meromorphic process X, let us consider 
the function 

tc(x) = 1 {X<0} £ ajePF, (5.36) 

where the coefficients aj and pj are positive and pj increase to +°° as j — > +°°. As 
was shown in [64], convergence of the series 

7>1 Pj 

implies convergence of the integral 

/ x 2 7t(x)dx < oo, 

./ — oo 

thus n(x) can be used to define the density of the Levy measure. Note that by con- 
stuction 7r(— x) is a completely monotone function. 

Using the Levy-Khintchine formula we find that the Laplace exponent is given 

by 

Y(z)= l -o 2 z 2 + Hz+z 2 Z " J ' , zeC, (5.37) 

and we see that y/(z) is a meromorpic function which has only negative poles at 
points z = — Pj- From the general theory we know that for q > equation \j/(z) = q 
has a unique solution z = <P(q) in the half -plane Re(z) > 0, and we also know that 
this solution is real. It can be proven (see [64]) that the same is true for equation 
z) = q- For q > all the solutions z — C, of y/(— z) =qin the halfplane Re(z) > 
are real and they satisfy the interlacing property 
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0<Ci<Pi<C2<p2<-... (5.38) 

When q = and y/(0) = E[Xi] < we have £i = 0, otherwise £i > 0. 

The following proposition gives an explicit formula for the potential density, gen- 
eralizing (5.35), which is expressed in terms of the roots ^ and the first derivative 
of the Laplace exponent. 

Proposition 5.5. 

( i) IfX is Meromorphic, then for all q > function ( x ) f s an infinite mixture 
of exponential functions with positive coefficients 

fi {q) (x) = - 1 £- rT - x>0. (5.39) 
j=\ r \ hj) 

(ii) If for some q > function u^\x) is an infinite mixture of exponential func- 
tions with positive coefficients, then X is a Meromorphic process. 

The first statement of the above proposition can be proved using Proposition 
5.1 and standard analytical techniques. See for example Corollary 2 and Remark 
2 in [65] as well as [66]. The second statement can be established using the same 
technique as in the proof of Theorem 1 in [64]. In both cases we leave all the details 
to the reader. 

Proposition 5.5 shows that we can compute the potential density and the scale 
function very easily, provided that we know and There seems little 

hope to find explicitly for a general q > and hence these quantities have to be 
computed numerically, which in turn requires multiple evaluations of y/(z). Com- 
puting \jf(z) with the help of the partial fraction decomposition (5.37) is not the best 
way to do it, as in general the series will converge rather slowly. Therefore it is 
important to find examples of Meromorphic processes for which y/(z) can be com- 
puted explicitly. Below we present several such examples, the details can be found 
in [63, 64]. 

(i) 0-process with parameter X € {3/2,5/2} 

W ( z ) = l -o 2 z 2 + iiz + c(-l) A - 1 / 2 (a +z/py- J coth (n^Ja + z//j) 

- c(-l) x - 1/2 a l - l coth(7ty/a) . (5.40) 

(ii) j3-process with parameter A e (1 , 2) U (2, 3) 

VA(z) = ^(7 2 z 2 + Mz + cB(l + a + z//3,l-A)-cB(l + a,l-A) , (5.41) 

where B(x,y) = r (x)T (y) / T (x + y) is the Beta function. 

The admissible set of parameters is a > 0, jU e R, c > 0, a > and j3 > 0. The 
parameters aj and Pj, which define the Levy measure via (5.36), are given as follows 
(see [64, 63]): in the case of 0-process we have 
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a J =lcpj 2X - 1 , Pj =p(a+j 2 ), (5.42) 
and in the case of j3 -process 

<tj = cp( J+ .^~ 2 y pj = P(a + j). (5.43) 

In the case of j3 -process we also have an explicit formula for the density of the Levy 
measure, 

e (l+a)px 

Moreover, it can be shown that the density of the Levy measure n{x) satisfies (up 
to a multiplicative constant) 

%{x) ~ \x\~ X , as x -> 0" , 
n{x) ~ S l+a )\ as jc->-oo. 

In particular, we see that the Levy measure always has exponential tails and the rate 
of decay of n(x) is controlled by parameters a and p\ See [64, 63]. The parameter 
c controls the overall "intensity" of jumps, while X is responsible for the behavior 
of small jumps: if A G (1,2) (resp. X G (2,3)) then the jump part of the process is 
of infinite activity and finite (resp. infinite) variation. 

The beta family of processes is more general than the theta family, as it allows for 
a greater range of the parameter X, which gives us more flexibility in modeling the 
behavior of small jumps. However, processes in the theta family have the advantage 
that the Laplace exponent is given in terms of elementary functions, which helps to 
implement faster numerical algorithms. 

Computing the scale function for Meromorphic processes is very simple. The 
first step is to compute the values of which are defined as the solutions to 
z) = q. The interlacing property (5.38) gives us left/right bounds for each £y 
(recall that Pj are known explicitly), thus we know that on each interval (p 7 ,p,+i) 
there is a unique solution to z) = q, and this solution can be found very effi- 
ciently using bisection/Newton's method. 

There is one additional trick that can greatly reduce the computation time for this 
first step. We will explain the main idea on the example of a j3 -process. Formula 
(5.43) shows that the spacing between the poles pj is constant and is equal to p\ This 
fact and the behavior of y/(z) as Re(z) — > °° implies that for j large, the difference 
— £y would also be very close to p\ Therefore, once we have computed for 
a sufficiently large j, we can use £y + j3 as a starting point for our search for 
In practice, starting Newton's method from this point gives the required accuracy 
in just one or two iterations. A corresponding strategy can be developed for the 9- 
processes: one should use the fact that y/pj — a/3 have constant spacing equal to 

VP- 
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Once we have precomputed the values of we compute the coefficients 
which is an easy task since we have an explicit formula for (z) . Now we can eval- 
uate u^ q \x) by truncating the infinite series in (5.39). Note that the infinite series 
converges exponentially fast (and much faster in the case of 9 -process), so unless x 
is very small, one really needs just a few terms to have good precision. On the other 
hand, when x is extremely small, it is better to compute the scale function by using 
the information about W^(0 + ) and W (<?) '(0+) given in Lemmas 3.1 and 3.2 and 
applying interpolation. 



5.6 Numerical examples 

In this section we present the results of several numerical experiments. Before we 
start discussing the details of these experiments, let us describe the computing en- 
vironment. All the code was written in Fortran90 and compiled with the help of 
Intel® Fortran compiler. For the methods which require multi-precision arithmetic 
we have the following two options. The first one is the quad data type (standard on 
Fortran90), which uses 128 bits to represent a real number and allows for the preci- 
sion of approximately 32 decimal digits. Recall that the double data type (standard 
on C/C++/Fortran/Matlab and other programming languages) uses 64 bits and gives 
approximately 16 decimal digits. The second option is to use MPFUN90 multi- 
precision library developed by Bailey [10]. This is an excellent library for Fortran90, 
which is very efficient and easy to use. It allows one to perform computations with 
precision of thousands of digits, though in our experiments we've never used more 
than two hundred digits. All the computations were performed on a standard 2008 
laptop (Intel Core 2 Duo 2.5 GHz processor and 3 GB of RAM). 

In our first numerical experiment we will compare the performance of the four 
methods described in Sections 5.2 and 5.3 in the case of a 0-process with X = 3/2 
(process with jumps of finite variation and infinite activity), whose Laplace exponent 
is given by (5.40). We will consider two cases a = or a = 0.25 and will fix the 
other parameters as follows 

ju = 2, c = 1, a = 1, j8 =0.5. 

Moreover, everywhere in this section we fix q = 0.5. 

As the benchmark for this experiment, we compute the values of the scale func- 
tion W^\x) for one hundred values of x, given by x ; - = //20, i— 1,2,..., 100. This 
part is done using the algorithm described in Section 5.6. The infinite series (5.39) is 
truncated at j = 1000 and the system precision is set at 200 digits (using MPFUN90 
library). This guarantees that our benchmark is within 1.0e-150 of the exact value. 

Next we compute the approximation (x) using Filon/ Gaver-Stehfest/ Euler/ 
Talbot algorithms for the same set of values of x. As the measure of accuracy of the 
approximation we will take the maximum of the relative error 
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relative error = max 

i<;<ioo 



\wM(xi)-WW(xi)\ 
W(i)( Xi ) 



(5.44) 



Let us specify the parameters for the Filon's method, as described in Section 5.2. 
We truncate the integral at b = 10 7 (or b = 10 9 ) in (5.25) and divide the domain 
of integration into ten sub-intervals, i.e. n — 10. Moreover, we set bj — b(j/n) 5 for 
j = 1,2, ... ,n so that there are more discretization points close to u = 0. The number 
of discretization points per each sub-interval [bj,bj+\] is fixed atN x = 10 4 (or N x = 
10 5 ). Note that Filon's method as described in Section 5.2 produces N x values of 
spacing between the points in x-domain is chosen at 8 X = 5/N x . In this 
way the grid x m = m8 x covers the whole interval [0,5]. In order to perform the Fast 
Fourier Transform we use Intel® MKL library. 



algorithm 


rel. error 


time (seconds) 


N x 


M 


system precision 


Filon (b = 10 7 ) 


3.0e-12 


0.18 


10 000 




double (64 bits) 


Gaver-Stehfest 


9.2e-12 


0.4 


100 


20 


44 (MPFUN90) 


Euler 


9.0e-14 


0.05 


100 


20 


quad (128 bits) 


Talbot 


2.1e-13 


0.025 


100 


20 


quad (128 bits) 


Gaver-Stehfest 


2.9e-21 


1.5 


100 


40 


88 (MPFUN90) 


Euler 


2.7e-25 


1.5 


100 


40 


40 (MPFUN90) 


Talbot 


2.4e-25 


0.7 


100 


40 


40 (MPFUN90) 


Gaver-Stehfest 


2.2e-39 


6.3 


100 


80 


176 (MPFUN90) 


Euler 


2.1e-48 


5.0 


100 


80 


80 (MPFUN90) 


Talbot 


3.1e-49 


2.5 


100 


80 


80 (MPFUN90) 



Table 1: Computing W (e >\x) for the 9 -process with a = 0.25. 



The results of this numerical experiment for computing W^\x) are presented 
in Table 1 . We see that each of the Gaver-Stehfest/Euler/Talbot algorithms produce 
very accurate results. When M = 40 we obtain more than twenty significant dig- 
its and when M = 80 we get almost fifty significant digits. The computation time 
varies with each algorithm and also depends on parameters. Note that when we use 
the quad date type for Euler/Talbot methods with M = 20, the computations are very 
fast, on the order of one hundredth of a second. When we use MPFUN90 library the 
computation time increases significantly. This is due to the fact that quad is a native 
data type in Fortran90 and computations done with an external multi-precision li- 
brary are slower. Next we see that Talbot algorithm is always about two times faster 
than Euler method. This is due to the fact that the finite sum in (5.29) has M terms 
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while Euler method (5.28) requires summation of 2M + 1 terms. We would like to 
stress again that the time values presented in this table are for computing the scale 
function for N x different values of x. For example, if one wants to compute the scale 
function for a single value of x using Talbot method with M = 40, it would take just 
seven milliseconds. 

Filon's method also performs well, computing the scale function for 10000 values 
of x with accuracy of around 12 decimal digits in just 0.18 seconds. As we see from 
the table, the Talbot algorithm appears to be the most efficient. This is not surprising 
given the discussion in Section 5.3.3 where it was pointed out that this algorithm 
is well suited for processes with completely monotone jumps (see also results by 
Albrecher, Avram and Kortschak [5] on ruin probabilities for completely monotone 
claim distributions). 



algorithm 


rel. error 


time (seconds) 


N x 


M 


system precision 


Filon (b = 10 9 ) 


8.0e-5 


2.9 


10 s 




double (64 bits) 


Gaver-Stehfest 


4.6e-12 


0.45 


100 


20 


44 (MPFUN90) 


Euler 


1.7e-13 


0.05 


100 


20 


quad (128 bits) 


Talbot 


1.5e-12 


0.023 


100 


20 


quad (128 bits) 


Gaver-Stehfest 


7.6e-21 


1.5 


100 


40 


88 (MPFUN90) 


Euler 


4.3e-25 


1.4 


100 


40 


40 (MPFUN90) 


Talbot 


2.9e-24 


0.71 


100 


40 


40 (MPFUN90) 


Gaver-Stehfest 


1.7e-38 


6.3 


100 


80 


176 (MPFUN90) 


Euler 


2.8e-48 


4.9 


100 


80 


80 (MPFUN90) 


Talbot 


9.1e-48 


2.4 


100 


80 


80 (MPFUN90) 



Table 2: Computing W^'(x) for the 0-process with a = 0. 



In Table 2 we present the results of computing the first derivative of the scale 
function. The results are very similar to those presented in Table 1, the major differ- 
ence now is that Filon's method is not producing a very good accuracy. This happens 
because the integrand in (5.15) decays very slowly (see Proposition 5.3), and even 
though we have increased b and N x , we still get only five significant digits. 

For our second numerical experiment we will consider a process with jumps of 
rational transform, which are not completely monotone. We define the density of 
the Levy measure as follows 



7r W = 1 {x<o} e *(l+cos(ax)). 



(5.45) 
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Note that this is an example of a process which has jumps of rational transform but 
whose distribution does not belong to the phase-type class. Recall that any distri- 
bution from the phase-type class can be characterized as the time to absorption in a 
finite state Markov chain with one absorbing state when the chain is started from a 
particular state. One can easily see that 7c(—x) can not be the density of such an ab- 
sorption time, since it is equal to zero forx = (2k + \)%/(a). The Laplace exponent 
of this process can be computed using formula (5.31) 

/ n 1 2 2 1 Z+ 1 , 1 

V(z) = «<7V + MzH r + 7 ttt 9 



2 ^ z+1 (z+l) 2 + a 2 1+a 2 ' 

We fix the parameters at 

d = 0.25, Ll = 2, a=A. 

For this set of parameters we have <P(q) ~ 0.37519 (recall that we have fixed q = 0.5 
throughout this section) and all the roots of yr(— z) = q are simple and are located at 

Ci, 2 ~ 0.97265 ± 4.05 18i, £ 3 ~ 0.64448, £ 4 ~ 64.7854 



algorithm 


rel. error 


time (seconds) 


N x 


M 


system precision 


Filon (b = 10 7 ) 


1.2e-ll 


0.14 


10 000 




double (64 bits) 


Gaver-Stehfest 


1.2e-5 


0.06 


100 


20 


44 (MPFUN90) 


Euler 


9.1e-14 


0.048 


100 


20 


quad (128 bits) 


Talbot 


9.2e-5 


0.022 


100 


20 


quad (128 bits) 


Gaver-Stehfest 


5.5e-10 


0.23 


100 


40 


88 (MPFUN90) 


Euler 


2.6e-25 


0.44 


100 


40 


40 (MPFUN90) 


Talbot 


2.8e-13 


0.22 


100 


40 


40 (MPFUN90) 


Gaver-Stehfest 


1.5e-27 


1.2 


100 


80 


176 (MPFUN90) 


Euler 


2.0e-48 


1.2 


100 


80 


80 (MPFUN90) 


Talbot 


2.9e-49 


0.6 


100 


80 


80 (MPFUN90) 



Table 3: Computing (x) for the process with jumps of rational transform. 



Once we have these roots the potential density u^ q \x) can be computed from 
(5.35) and the scale function using relation (5.4). The parameters for Filon's method 
are the same as in our previous experiment with 0-process. The results are presented 
in Table 3 and are seen to be similar to those of the 0-process, with the exception 
that the Talbot method is not performing as well for small values of M. This is most 
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likely due to the fact that function (z) has non-real singularities in the half -plane 
Re(z) < and, as discussed in Section 5.3.3, this may cause problems for the Talbot 
method. However, when M — 80, the Talbot method is performing very well just as 
it did for the case of the 0-process. 

For our final numerical experiment, we consider a spectrally-negative Levy pro- 
cess with the Levy measure defined by 



n(dx) = cil {x<0} 



r dx + c 2 l { , <0} — -dx 



(5.46) 



C3 



8- a3 (dx) + c 4 l {x< _ aA} e x dx. 



Here the admissible range of parameters is c,- > 0, A; > 0, (X\ e (— °°,2) \ {0, 1}, 
(X2 € (0, 2) \ { 1 } and a, > 0. Using the Levy-Khintchine formula and standard results 
on tempered stable processes (see proposition 4.2 in [31]) we find that the Laplace 
exponent of X can be computed as follows 



1 



y(z) = ^o 2 z 2 + nz + c 1 r(-a l ) [(z + Xi) ai - A" 1 -za x X[ 



(5.47) 



+ c 2 r(-a 2 )z a2 + c 3 [e-" 3Z - 1] + c 4 



z+l 





a 




Cl 


h 


ai 




a 2 


C3 


a 3 


c 4 


«4 


Set 1 





2 


1 


1 


0.5 


1 


1.5 












Set 2 





2 


0.05 


1 


-5 


0.25 


1.5 












Set 3 


0.25 


2 


1 


1 


0.5 












1 


1 


Set 4 





2 


1 


1 


0.5 







1 


1 








Table 4: Parameter sets. 



We see that the Levy measure (5.46) is a mixture of (a) the Levy measure of 
a tempered stable process, (b) the Levy measure of a stable process, (c) an atom 
at — «3 and (d) an exponential jump distribution shifted by —04. Table 4 specifies 
parameter choices for four different sets of Levy processes which capture the fol- 
lowing characteristics. 

Set 1 : Jumps of infinite variation with a completely monotone Levy density; no 
Gaussian component. 

Set 2: Jumps of infinite variation which have smooth, but not completely mono- 
tone Levy density; no Gaussian component. 

Set 3: Jumps of infinite activity, finite variation and discontinuous Levy density; 
non-zero Gaussian component. 
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Set 4: Jumps of infinite activity, finite variation, the Levy measure has an atom; 
no Gaussian component. 

In order to compare the performance of the algorithms we need to have a bench- 
mark. In this experiment we don't have any explicit formula for the scale function, 
thus we have to compute the benchmark numerically. We have chosen to use Filon's 
method for this purpose as it is quite robust. For each case we have tried different pa- 
rameters to ensure that we have at least 10-11 digits of accuracy. For example, for the 
parameter Set 4 we fix b — 10 8 , c = 0.1, divide the interval of integration in (5.25) 
into 1000 subintervals and set bj ~ b(j/n) 2 for j — 1,2, ... ,n and n = 1000. Finally 
we fix the number of discretization points for each subinterval at N = N x = 5 x 10 6 . 

As in our previous numerical experiment, we compute the scale function for 
one hundred values of x, given by x, = z'/20, i — 1,2, ... , 100. The relative errors 
presented below are maximum relative errors over this range of x-values. 



algorithm 


rel. error 


time (seconds) 


N x 


M 


system precision 


Filon (b = 10 7 ) 


2.6e-10 


0.46 


10 000 




double (64 bits) 


Gaver-Stehfest 


1.2e-12 


1.1 


100 


20 


44 (MPFUN90) 


Euler 


1.3e-12 


0.05 


100 


20 


quad (128 bits) 


Talbot 


7.3e-13 


0.024 


100 


20 


quad (128 bits) 



Table 5: Computing W^(x) for the process with parameter Set 1. 



First we compute the scale function for the parameter Set 1 and present the results 
in Table 5. In this case we have a process with completely monotone jumps, and we 
see that the situation is exactly the same as in our earlier experiment with 0-process. 
All four methods provide good accuracy, although the Talbot algorithm is faster and 
gives better accuracy. 



algorithm 


rel. error 


time (seconds) 


N x 


M 


system precision 


Filon (b = 10 7 ) 


8.4e-10 


0.47 


10 000 




double (64 bits) 


Gaver-Stehfest 


3.9e-12 


1.1 


100 


20 


44 (MPFUN90) 


Euler 


3.9e-12 


0.05 


100 


20 


quad (128 bits) 


Talbot 


3.6e-12 


0.024 


100 


20 


quad (128 bits) 



Table 6: Computing (x) for the process with parameter Set 2. 
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In Table 6 we show the results for the parameter Set 2. In this case the jumps 
have smooth, but not completely monotone density. However, this does not seem to 
affect the results, and the performance of all four algorithms is similar to the case of 
parameter Set 1 . 



algorithm 


rel. error 


time (seconds) 


N x 


M 


system precision 


Filon (b = 10 7 ) 


4.6e-ll 


0.46 


10 000 




double (64 bits) 


Gaver-Stehfest 


3.3e-5 


1.1 


100 


20 


44 (MPFUN90) 


Euler 


1.4e-6 


0.06 


100 


20 


quad (128 bits) 


Talbot 


5.0e-3 


0.028 


100 


20 


quad (128 bits) 


Gaver-Stehfest 


7.0e-6 


4.5 


100 


40 


88 (MPFUN90) 


Euler 


9.1e-7 


4.7 


100 


40 


40 (MPFUN90) 


Talbot 


1.7e-3 


2.3 


100 


40 


40 (MPFUN90) 


Gaver-Stehfest 


l.le-6 


22 


100 


80 


176 (MPFUN90) 


Euler 


3.0e-6 


18 


100 


80 


80 (MPFUN90) 


Talbot 


2.9e-4 


8.8 


100 


80 


80 (MPFUN90) 



Table 7: Computing W^ q \x) for the process with parameter Set 3. 



In Table 7 we see the first example where the jump density (and therefore the 
scale function) is non-smooth. In this case we have W^(x) e C 3 (0,°°) (see The- 
orem 3.11), and we would expect that the discontinuity in the fourth derivative of 
W( q \x) shouldn't affect the performance of our numerical methods. Unfortunately, 
this is not the case. We see that, while Filon's method is still performing well, the 
Talbot algorithm provides only 2-3 digits of precision. The Gaver-Stehfest and Eu- 
ler algorithms do a slightly better job and produce 5-6 digits of precision. However, 
when we increase M we do not see a significant decrease in the error, and, in fact, 
it is not clear whether these three methods would converge to the right values at all. 
This behaviour, specifically problems with Laplace inversion when the target func- 
tion is non-smooth, is quite well-known. See for example the discussion in [1] and 
Section 14 of [3]. 

Finally, in Table 8 we present the most extreme case of parameter Set 4. In this 
case the scale function is not even continuously differentiable, so we should expect 
to see some interesting behavior of our numerical algorithms. We see that Filon's 
method is still performing very well. Indeed, the maximum relative error is around 
5.3e-8. The error achieves its maximum at x = 1, which is not surprising since this 
is the point where W^ q \x) is not differentiable (see figure 1). If we remove this 
point, then the maximum relative error drops down to 4.6e- 1 1 , which is an extremely 
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algorithm 


rel. error 


time (seconds) 


N x 


M 


system precision 


Filon (b = 10 7 ) 


5.3e-8 (4.6e-ll) 


5.1 


10 s 




double (64 bits) 


Gaver-Stehfest 


1.6e-3 


1.1 


100 


20 


44 (MPFUN90) 


Euler 


3.2e-4 


0.05 


100 


20 


quad (128 bits) 


Talbot 


4.6e-3 


0.024 


100 


20 


quad (128 bits) 


Gaver-Stehfest 


8.0e-4 


4.5 


100 


40 


88 (MPFUN90) 


Euler 


3.2e-4 


4.6 


100 


40 


40 (MPFUN90) 


Talbot 


2.4e-3 


2.2 


100 


40 


40 (MPFUN90) 


Gaver-Stehfest 


3.9e-4 


22 


100 


80 


176 (MPFUN90) 


Euler 


8.1e-3 


18 


100 


80 


80 (MPFUN90) 


Talbot 


1.2e-3 


8.7 


100 


80 


80 (MPFUN90) 



Table 8: Computing (x) for the process with parameter Set 4. 



good degree of accuracy for such non-smooth function. By contrast, the other three 
algorithms are all struggling to produce even four digits of precision. Again, it is not 
clear whether these algorithms converge as M increases. In particular we see that the 
error of the Euler method appears to increase. 

We obtain an overview of what is happening in this case by inspecting the graph 
of W^'(x) (see Figure 1). We see that W^'(x) has a jump at x = 1, which is the 
point where the Levy measure has an atom. This behavior is expected due to Corol- 
lary 2.5. Then we observe that while W^'(x) is continuous at all points 1, the 
second derivative W^"(x) has a jump at x = 2. Although there is no theoretical 
justification for the appearance of such a jump in this text, this observation is con- 
sistent with the results on discontinuities in higher derivatives of potential measures 
of subordinators obtained by Doring and Savov [35]. 



5.7 Conclusion 

In the previous section we have presented results of several numerical experiments, 
which allow us to compare the performance of four different methods for computing 
the scale function. 

Filon's method has the advantage of being the most robust. We are guaranteed to 
have an approximation which is within 0(N~ 3 ) of the exact value, provided that we 
have chosen the cutoff b to be large enough. Another advantage of this method is that 
it gives not just a single value of W^ q \x), but an array of N values W^ q \xj) (where 
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(a) fW'(j) (b) log 10 (relative error) with M= 20 



Fig. 1: Computing W^ q ' (x) for the process with parameter Set 4. The relative errors 
produced by the Filon, Gaver-Stehfest, Euler and Talbot methods correspond to the 
green, black, blue and red curves respectively. 



the xi are equally spaced) at the computational cost of just 0(Nln(N)) operations. 
This means that by increasing N we get a better accuracy and, at the same time, 
more values of W^ q '(x), and the computational time would increase almost linearly 
in N. This is different from the other three methods, where the computational time 
needed to produce N values of (x) is 0(NM). 

At the same time, there are several disadvantages relating to Filon's method. 
First of all, this algorithm requires quite some effort to program. Secondly, when 
the integrand decreases very slowly, it may be necessary to take an extremely large 
value of the cutoff value b. This would have to be compensated by an increase in 
the number of discretization points N (see Table 2). Finally, unlike the other three 
methods, which depend on a single parameter M, Filon's method has two important 
parameters, b and N (as well as bj and n if one uses (5.25)), and it is not a priori 
clear how to choose these parameters. Usually, reasonable values of b and N can be 
found after some experimentation, which is of course time consuming. 

The essential property of the other three methods, i.e the Gaver-Stehfest, Euler 
and Talbot algorithms, is that all of them require multi-precision arithmetic. If the 
scale function is smooth and we need just 10-12 digits of precision, then we can use 
the quad (128 bits) precision, which is available in some programming languages, 
such as Fortran90. In all other cases one has to use specialized software which is ca- 
pable of working with high-precision numbers. Whilst very efficient multi-precision 
arithmetic libraries do exist, such as MPFUN90 by D.H. Bailey [10], they are usu- 
ally much slower than the native double/quad precision of any programming lan- 
guage. 

An advantage of the Gaver-Stehfest/Euler/Talbot methods is their simplicity. 
These algorithms are very easy to program, and they all depend only on just a single 
integer parameter M. 

The performance of the the Gaver-Stehfest, Euler and Talbot algorithms depends 
a lot on the smoothness of the scale function W^ q \x). When W^ q \x) is not suf- 
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ficiently smooth it might be better to use Filon's method. When the process has 
completely monotone jumps, Talbot's algorithm is arguably preferable to the other 
methods as it produces excellent accuracy. Moreover, it is almost twice as fast as 
its nearest competitor, the Euler algorithm. When the process has a smooth, but not 
completely monotone, jump density the Euler algorithm is recommendable. One 
can still use the Talbot algorithm in this case, however, this should be done with a 
caution as it is not guaranteed to work (see the discussion in Section 5.3.3 and the 
results of our second numerical experiment). 

Finally, the Gaver-Stehfest algorithm is somewhat special. It is the only algorithm 
which does not require the evaluation of \jf(z) at complex values of z. Our results 
show that this method is not the fastest and not the most precise, but it can still 
be useful in some situations. For example, this method might be the only available 
option when the Laplace exponent \jf(z) is given by a very complicated formula 
which is valid for real values of z and, moreover, it is not clear how to analytically 
continue it into the complex half-plane Re(z) > 0. The Gaver-Stehfest algorithm can 
also be helpful when one is using multi-precision software which does not support 
operations on complex numbers. 
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